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(ABSTRACT) 

The Thermal Radiation Group at Virginia Polytechnic Institute and State University has 
been working closely with scientists and engineers at NASA’s Langley Research Center 
to develop accurate analytical and numerical models suitable for designing next- 
generation thin-film thermal radiation detectors for earth radiation budget measurement 
applications. The current study provides an analytical model of the notional thermal 
radiation detector that takes into account thermal transport phenomena, such as the 
contact resistance between the layers of the detector, and is suitable for use in parameter 
estimation. It was found that the responsivity of the detector can increase significantly 
due to the presence of contact resistance between the layers of the detector. Also 
presented is the effect of doping the thermal impedance layer of the detector with 
conducting particles in order to electrically link the two junctions of the detector. It was 
found that the responsivity and the time response of the doped detector decrease 
significantly in this case. The corresponding decrease of the electrical resistance of the 
doped thermal impedance layer is not sufficient to significantly improve the electrical 
performance of the detector. Finally, the “roughness effect” is shown to be unable to 
explain the decrease in the thermal conductivity often reported for thin-film layers. 
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Chapter 1: Introduction 


For many years the Thermal Radiation Group (TRG) at Virginia Polytechnic Institute and 
State University, under the direction of Professor J.R. Malian, has been involved in 
developing advanced numerical tools for modeling the optical, thermal and electrical 
processes in thermal radiation detectors. During the last two years the Group, working 
with scientists and engineers at NASA’s Langley Research Center, has pursued the 
design of a new thermal radiation detector to be used for more accurate measurements of 
the earth radiation budget (ERB) from space. The ERB refers to the balance at the top of 
the atmosphere between the incoming energy from the sun and the outgoing thermal and 
reflected energy from the earth. 

The growing interest in ground-based measurements using instruments such as 
pyranometers [Smith, 1999], which are cheaper and easier to monitor than space-borne 
instruments, makes it crucial for the proponents of space-borne measurements to develop 
more reliable instruments using the latest technologies available in order to remain 
competitive. The sensors of these latter instruments have to be as small as possible, with 
extremely short tune response and a maximum responsivity, so that they can be organized 
into linear or focal-plane arrays [Weckmann, 1997; Barreto, 1998; Sanchez, 1998; 
Sorensen, 1998]. These factors all favor the development of a thin- film multilayer 
thermal radiation detector technology, as illustrated in Figure 1. In this approach six 
layers of material are sputtered onto a substrate to form the detector. This allows layers 
on the order of a few microns thick. The absorption of thermal radiation by the absorber 
layer increases the temperature of the active junction, and the difference of temperature 
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between this active junction and the reference junction generates an electromotive force 
proportional to the temperature difference between them. 

Analysis of heat transfer in the resulting thin-film thermal radiation detector is essential 
for its successful design and for the prediction of its performance. Creation of the heat 
transfer model requires applying either macroscale or microscale transport theories, 
depending on the physical dimensions (length scales) of the detector and its frequency 
response characteristics (time scale). If macroscale heat transfer theory is applied to a 
microscale problem in a situation for which it is inappropriate, then a significant error in 
the calculated heat transfer rate or temperature distribution can result [Flik et al., 1992], 
The model also requires the precise knowledge of the thermophysical properties of 
materials used in the conception of the detector. With the contemplated geometry and 
fabrication technique, several interesting questions arise, and the performance of the 
detector will depend on the response to these questions. Ignoring these nuances could 
lead to very poor thermal models of the detector. 

Normally the thermal resistance offered to heat conduction by a uniform layer is directly 
proportional to its thickness. When two materials having different thermal conductivities 
are in mechanical contact, a temperature discontinuity can occur at their interface. Tins 
temperature jump results from contact resistance between the two layers. In fact, even 
with thin-films, perfect contact at the interface occurs only at a limited number of spots 
[Hmina et al., 1997; Kelkar et al., 1996] and the void found elsewhere between the layers 
is filled with gas or a vacuum. It is possible that the intrinsic thermal resistance of a thin 
film is significantly smaller than the contact resistance at the interface. Modeling the 
detector without taking into account this hypothetical contact resistance could lead to a 
serious departure from reality. Even if the values of the contact resistances between the 
different layers are not currently known, it is important to include them in the thermal 
models. These more general models can then be used by other members of the Thermal 
Radiation Group who are developing parameter estimation techniques aimed at 
recovering the thermophysical properties of the detector. The use of fictitious contact 
resistance values in the model will reveal the importance of that factor to the performance 
of the detector. 
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The correct treatment of the thin-film effect is another key factor in the process of 
obtaining a good thermal model. The thin -film effect has been defined in terms of the 
departure of the thermal conductivity of a layer from its bulk value as the thickness of the 
layer approaches the order of magnitude of the mean free path of the energy carriers, a 
definition that presupposes that Fourier’s law of heat conduction governs thin-film heat 
transfer. Flowever the thin- film effect implies the use of microscale heat transfer theories. 
In these theories the classical Fourier’s law no longer describes the temperature and heat 
flow fields, and the thennal conductivity loses its meaning as an intrinsic thermal 
property of the thin-film material. Several analytical [Flik et al., 1990; and Kumar et al., 
1994] and experimental [Nath et ah, 1974; Lambropoulos et ah, 1989; Lee et ah, 1997; 
and Oram et al., 1998] studies have shown that the thennal conductivity of thin-film 
layers, defined as the heat flux divided by the temperature difference across the layer per 
unit thickness, can be several orders of magnitude lower than that of the bulk material. 
The closed- fonn models for computing thennal transport in such cases require the 
knowledge of material properties such as the mean fine path that are unknown for 
materials contemplated for use in the notional detector. 

In the cunent effort the use of a geometrical property characterizing the mechanical 
roughness of the film layer is considered in an attempt to show the dependence between 
the effective thennal conductivity and the layer thickness. It is evident that in a thin film, 
the roughness of the layer can be of the same order of magnitude as the layer thickness 
itself. The idea here is to try to explain the decrease of the thermal conductivity often 
attributed to the thin-film effect in tenns of a roughness effect. A numerical thennal 
model derived from a macroscale approach with the roughness as a parameter is used to 
analyze the dependence of the effective thennal conductivity on the roughness of the 
layer. 

The overall distance between the two thennocouple junctions of the notional detector 
shown in Figure 1 is less than fifty microns. Another problem that arises is how to 
electrically connect these two junctions. Fabrication and size limitations preclude the use 
of a traditional electrical conductor (a “wire” lead) to connect these two junctions. The 
NASA design team led by Kist [1999] has proposed that the thermal impedance layer 
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between the two thermocouple junctions also double as an electrical conductor. One idea 
would be to dope the thermal impedance layer with carbon particles. It is hoped that these 
particles would then effect the electrical connection between the junctions without 
completely eliminating the thermal impedance needed to maintain the temperature 
difference between the active and reference thermocouple junctions. Two different 
approaches are considered in order to characterize the thermal and electrical behavior of 
this new detector configuration. The first approach consists of using a model adapted 
from Vick and Scott [1998]. The second approach models the heterogeneous thermal 
impedance layer of the detector using micromechanical approximations to obtain 
averaged electrical and thermophysical properties. In the first approach the thermal 
impedance layer is modeled as a matrix with embedded small particles, and heat is 
assumed to be transferred from the matrix to the separated particles. The temperature 
distribution within such a layer is obtained by solving two coupled equations derived 
from the energy balance of the matrix and the particles. In the second approach, an 
equivalent thermal conductivity and heat capacity are used to represent an appropriately 
layered medium, and the detector is modeled as if the layers were homogeneous. The 
same homogeneous approximations are also used in order to analyze the change of the 
electrical resistance of the doped thermal impedance layer. In either case, the resulting 
model would need to be suitable for later use in a parameter estimation scheme. 

Objectives 

The objectives of the current investigation are to explore the following three questions: 

• How might the presence of contact resistance influence the response of the 
detector? 

• Can the so-called thin-film effect be explained and treated as a roughness effect? 

• What is the effect of the addition of electrically conducting particles in the 
thermal impedance layer on the performance of the detector? 

Even if the answer to the second question is negative, a secondary goal is to advance the 
Group’s knowledge of this interesting phenomenon. 
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In pursuing the answers of these questions the author has profited from the fruits of 
earlier investigations by members of the Thermal Radiation Group. Weckmann’s [1997] 
model of the dynamic electrothermal response of a similar thermal radiation detector 
ignored contact resistance and did not speculate on the possible consequences of a “thin- 
fUm” effect. Sanchez [1998] formulated a Monte-Carlo ray-trace radiative model of an 
integrated detector concept in which a linear array of the detectors from Weckmann’s 
thesis was bonded to the wall of a V-groove cavity, and Barreto [1998] designed 
experiments for estimating the electro-thermophysical properties of the materials in 
Weckmann’s detector concept. More recently Sorensen [1998] performed analytical and 
experimental characterizations of some aspects of the detector. Sorensen is currently 
developing genetic and hybrid algorithms and designing experiments to estimate the 
electro-thermophysical properties of detectors such as the one shown in Figure 1 . 

Chapter 2 addresses the problem of modeling the effects of contact resistance for use in 
the optimal design of experiments aimed at recovering its value. 
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Chapter 2: The Effect Contact Resistance 


Three questions bearing on the electrothermal performance of thin-film thermal radiation 
detectors were posed in Chapter 1. In the current chapter, one of these questions, that 
concerning contact resistance, is explored in detail. 

2.1 Description of the notional thermal radiation detector 

The notional thermal radiation detector considered in this thesis consists of six parallel 
layers of different materials, as illustrated in Figure 1 . A 1 .0-pm thick platinum layer is 
deposited on an aluminum-nitride substrate, and then a l.O-pm zinc-antimonide layer is 
sputtered over the platinum. Next, a 25-pin thick thermal impedance layer is deposited 
over the zinc-antimonide layer. Then additional layers of zinc-antimonide and platinum 
are deposited on the thermal impedance layer. Finally, the platinum is coated with a 10- 
prn thick black absorber layer to increase the absorption of thermal radiation. The 
platinum/zinc- antimonide layer near the substrate forms the reference junction of the 
thermocouple, and the platinum/zinc-antimonide layer on the thermal impedance layer 
forms the active junction. The junction materials are chosen to provide the highest 
available Seebeck coefficient [Weckmann, 1997]. In this baseline configuration the 
thermal impedance is not doped with carbon particles. The doped configuration is 
considered in Chapter 4. The thermal impedance layer considered in the current chapter is 
assumed to be fabricated from Parylene. The nominal properties of the material used in 
the detector are given in Table 1. 
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Incident thermal radiation 



Absorber layer (10 (tm) 

Zinc- ant imonide layer (1 um) 

Platinum layer (1 jam) 

’Thermal impedance layer 
(25 


[Not to scale] 


Figure 1: Configuration of the notional thermal radiation detector considered in tins 
thesis. 
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Table 1: Thermophysical properties of materials used in the baseline thermal radiation 
detector model [Weckmann, 1997]. 



Mass 

density 

(kg/m 3 ) 

Specific 

heat 

( J/kg-K ) 

Conductivity 

(' W/mK ) 

Diffusivity 

(m 2 /s) 

Seebeck 

Coefficient 

(V/K) 

Aluminum 

Substrate 

3260 

800 

165 

6.33X10" 5 

N/A 

Absorber 

Layer 

1400 

669 

0.209 

2.23 xl0“ 7 

N/A 

Parylene 

1289 

712 

0.084 

9.15xl0~ 8 

N/A 

Platinum 

21450 

133 

71.6 

2.51 x 10" 5 

N/A 

Zinc 

Antimonide 

6880 

200 

60 

4.4x1 O' 5 

N/A 

Carbon 

2620 

710 

1.59 

8.55xl0~ 7 

N/A 

Platinum/Zinc- 
Ant imonide 
Junction 

N/A 

N/A 

65.3 

3.09X10" 5 

960x1 0" 6 
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2.2 One-dimensional model of the detector with contact resistance 

Due to the thickness of the layers of the thermal detector, on the order of a few microns, 
the thermal resistance resulting from the imperfect contact between layers, illustrated in 
Figure 2(a), may play an important role in the overall thermal resistance. In fact, a careful 
search of the literature provided no insight into the nature of the interface and the 
consequent contact resistance at the interface between the sputtered materials in the 
notional detector. With contact resistance, the temperature profile through the detector 
would be discontinuous at the interfaces between layers. The first objective of the model 
is to validate the results obtained assuming perfect thermal contact. A preliminary study 
is also carried out to evaluate the effect of the thickness of a Parylene thermal impedance 
layer on detector time response and responsivity. 

In order to develop a one-dimensional analytical model we assume for both the steady- 
state and transient cases that: 

• the thermophysical properties (conductivity, specific heat, density) are constant 
and uniform, 

• thermal radiation is negligible, 

• Fourier’s law is valid, and 

• internal heat generation is absent. 
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Active 

junction 


Reference 

junction 


(a) 



(b) 

[Not to scale] 



8 Absorber layer H Zinc-Antinomide layer EH void 

■ Platinum layer ■ Thermal impedance layer 

Figure 2: (a) Model of the detector with contact resistance, and (b) boundary conditions. 
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2.2.1 Steady-state temperature distribution 


For each layer (i) the steady-state temperature rise 9 s j (x) above the substrate temperature 
is governed by 


dx 2 


= 0 


x,.<x<x i+1 , i = 1,2,.. .,6, (2-1) 


subject to the boundary conditions, indicated in Figure 2(b), 


Os, i = o 


at x = 0 , 


r do 


k. —r~ = h, M (6 S , -d sM ) 

dx 


k dO^ = k dd sM 


v ; , ,v /+l 

dx 


dx 


at the interfaces .r = x . , , 


/ = 1,2,.. .,5, 


and 


, d0 St 6 

^6 — — = q , 

dx 


at x = x 6 . 


(2-2) 

(2-3) 


(2-4) 


(2-5) 


The heat flux q is uniform through the different layers. 

The differential equations given by Equations 2-1 through 2-5 can be solved directly. 
For the first layer 


d 2 e sX 

dx 2 


= 0 


( 2 - 6 ) 


9 s l =0 , at x = 0 , 


(2-7) 


and 


q=k x 


dO s ,i 

dx 


at x = jq . 


( 2 - 8 ) 
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The temperature profile is then 


( 2 - 9 ) 


S,.M = T x 


For the second layer 

d 2 6 


s,2 


dx 2 


= 0 , 


' K ( ^ L = h U2 [0 s l (x 2 ) - 0 2 (x ., )] at a: = Xy 

dx 


and 


q = k 2 


dd 


s, 2 


dx 


at x = x 0 


Thus, d s2 (a) = q 


X 


1 

+ AC, * 

& 2 

k , Au 

v 1 2 > 

+ — 
k l2 _ 


More generally the steady- state temperature distribution for any layer is given by 


and 


e sl =±x 
’ K 


d s , = q 


v JzL 

f + 5 >~ 

K i 7=2 


1 1 


k J - 1 k J 


m=i - 1 1 

+ 1 


“ h 


i > 2. 


w=l 1 


( 2 - 10 ) 

( 2 - 11 ) 

( 2 - 12 ) 

( 2 - 13 ) 

( 2 - 14 ) 
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2.2.2 Transient temperature distribution 


The differential equation describing the evolution of the transient temperature rise 
6 ] (x, t ) above the substrate temperature for each layer is 


d 2 0 j (x,t) _ 1 dO t (x.,t) 
dx 2 a, d t 


x,<x.<x. i+1 , t > 0 , 


(2-15) 


where i = 1,2,..., 6 , subject to the boundary conditions 
<9, (0, t) = 0 at 


r 



= K + M-0 m ) 


V 


k: 


dx i+1 dx 


and 


at the interfaces 
i' = l,2,.. .,5 


x = , (2-16) 

x = x M , (2-17) 

t >0 , (2-18) 


dO 

k 6 — ± = q atx = x 6 , t>0 • (2-19) 

dx 

The initial condition is 0 ! (x,0) = 0 . (2-20) 

This problem can be solved by different techniques, including the finite difference 
method, by the finite element method and by finite integral transforms. However an 
analytical solution being more convenient for later parameter estimation studies, the 
orthogonal expansion technique [Ozisik, 1993] for a multi-layer medium with perfect 
contact between layers is adapted to the transient model described by Equations 2-15 
through 2-20. The problem involves a nonhomogeneous boundary condition at x = x 6 (at 

the top of the detector) and so cannot be directly solved using separation of variables. In 
order to transform this nonhomogeneous boundary condition into a homogeneous one, we 
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consider that 9 : (x,t) is obtained by the superposition of two problems, as suggested in 
Figure 3: 

• a steady-state problem with specified heat flux at x = x 6 : 9 si (x) 

• a transient problem with insulated boundary at x = x 6 : 9 tj (x, t) 


Thus 


9 i (x,t) = 9 si (x) + 9 ti (x,t) . 


( 2 - 21 ) 


t X 

4 q d9 t6 (^6 »0 


1 

r ^r \r 'r ^r \ 

r i 


dx 


d 1 9 j (x,t) 1 d9 j (x,t) 


dx “ 

+ 

d 2 9 ri (x,t) _ 1 d9 ri (x,t) 

dx 2 a, dt 

9,(x, 0) = 0 

dx 2 a t dx 

9 t i (x,0) = -9 s t (x) 


^(0,0 = 0 ^,i(0) = 0 9 rl (0,t) = 0 


Figure 3: Illustration of the superposition principle. 


The steady- state problem is described by an ordinary differential equation and has 
already been solved; the solution is given by Equation 2-14. 

The transient solution 9 tj (x,t) is governed by a homogeneous partial differential equation 
and is solved using the separation of variables technique. Let us assume that#, , (x,t) may 
be expressed 

9 ti (x,t) = X,(x)m . (2-22) 
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Substitution into Equation 2-15 yields 


3 2 (X,(x)r(r)) _ 1 3(X,(s)r(0) 

dx 2 a , dt 


(2-23) 


Thus 


a i d 2 X j _ 1 dT 
X . dx 2 ~ T dt 


(2-24) 


where the eigenvalue fd is the separation constant. The minus sign in front of fd 1 is 
required to capture the exponential nature of the variation with time. The change of 
variables transforms the partial differential equation into two ordinary differential 
equations, 


and 


d£m 

dt 


+AX=o 


(2-25) 


d X, m B 2 

} > m , > rn y 

, 2 ' ^ ;',/w 

dx a, 


0 


(2-26) 


The subscript m is used in Equations 2—25 and 2-26 to indicate that there are an infinite 
number of eigenvalues f3 rn . 

For a given eigenvalue the solution of Equation 2-25 is obtained directly as 

F m (t) = e~^ t . (2-27) 
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The complete solution is then of the form 


0=£c,x,.,(i)c-® , (2-28) 

m = 1 


where the coefficients C m are eventually obtained by exploiting the orthogonahty 
property. 


The ordinary differential equation to be solved to obtain the eigenfunctions X jm and the 
corresponding eigenvalues is 


d 2 *,, fa 


x,<x<x M ,i = 1, ... , 6 , (2-29) 


with boundary conditions 


*i»(0) = 0 


at 


x = 0 


(2-30) 


r , dx i, m , dX i+lm 
k, = k M 


dx 


dX,, 


dx 


at the interfaces 


k i ~r L = h u+1 (X Km -X i+Um ), i = 1,2,.. .,5 
dx 


X = X: 


(2-31) 


(2-32) 


and 




dx 


0 


at .r = a'. 


(2-33) 


The general solution of this eigenvalue problem is [Ozisik, 1993] 




x ,Jx) = A> sin 


P, 


X 


( 


+ B, _ cos 


Pn 


X 


(2-34) 
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Let us assume that //, = 


(2-35) 


Then Equation 2-34 becomes 


*,>(*) = 4>sin //,— +B i cos 77 ,.— 

X X. 


(2-36) 


The derivative of Equation 2-36 is 


dX, ft f x ) . f x 

—^ = ~r= A> c °s A £,>,sin 77,— 

dx ^a, I X, I I X,. 


(2-37) 


Substitution of Equations 2-35 and 2-37 into Equations 2-31 and 2-32 yields 


A+i.m s in — +B, + , w cos —77, 


= [ sin (7, ) + u+ 1 cos(77, )] + 5,. [cos(77, ) - sin(77, )] 


(2-38) 


A + um cos - B i+l m sin — 77, +1 = K l i+l [A i m cos(77, ) - sin(77, )] , (2-39) 
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where 


(2-40) 


H u+l = 


h,»4* 


K a i+ i 
k i+ 1 V a i 


(2-41) 


In matrix form Equations 2-38 and 2-39 can be written 


sin(x,77, +1 / x M ) cos(v/, +l /x M ) T A+u 

cos(a: jj m / x, +1 ) - sin(jf, jj. +1 / x i+1 )J A + i,„ 


sin(/7 , ) + H u+l 0 m cos(;;, ) cos(;;, ) - //, sin(;;, )¥ A. n 

K u + iCos(;;,) sin(;;,) £,> 


Solving Equation 2-42 for the unknown coefficients A i m and B j m yields 


where 


(2-42) 


X X 

A + I, m = u ,, m sin —11m +V,_cos —tj i+1 


i = 1,2,. ..,5 


(2-43) 


X. X- 

= U,, m cos —r/ M ~V i m sin —Tj i+l , i = l,2,...,5 


(2-44) 


U Km = A. m [sin(//, ) + p m H. i+l cos(;/, )j+ B [cos(^ ) - H. i+1 fi m sin(;/, )J (2-45) 
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Vf, m = K u + 1 U, > co <n, ) - B i m sin (Tj i )j 


(2-46) 


The boundary condition given by Equation 2-30 requires that B l m = 0 . Thus we can set 

A, m = i • 

All of the coefficients A j m and B im can now be computed using Equations 2-43 and 2- 
44. Introducing A l m = 1 into Equation 2-42 and expanding it for the six layers of the 
detector yields 


Ml 

Ml 

M\ 

0 

0 

0 

0 

M\ 

M\ 

Ml 

0 

0 

0 

0 

0 

Ml 

Ml 

Ml 

Ml 

0 

0 

0 


Ml 

Ml 

Ml 

0 

0 

0 

0 

0 

Ml 

Ml 

Ml 

Mj 

0 

0 

0 

Ml 

Ml 

Ml 

Ml 

0 

0 

0 

0 

0 

Ml 

Ml 

0 

0 

0 

0 

0 

Ml 

Ml 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 



,(2-47) 


where 

M\ = sin^ +H l 2 cos rj l 
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M 2 =K 12 cos tj x 



Ml = sin rj 2 + H 23 cos r / 2 
Ml = cos r / 2 - H 2 3 sin r / 2 



Ml = K 2 3 cos t ] 2 
M l — —K 2 3 sin 77 2 




M s 4 = sin + H 34 cos 
M 5 5 = cosr/ 3 -H 34 sin^ 3 
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M: = - sin 


^3 ' 

— Va 

^4 , 


M 5 = -cos 


f x. ' 






Ml =K i4 cosT] i 


M l =-^3 4 Sin7 73 


Ml = -cos 


r x 3 

—r/4 

K X 4 , 


M 6 = sin 


^3 ' 

— ^4 

y 


M* = sin/7 4 +// 45 


M] — COS Tj 4 H 4 5 


M 7 = - sin 


'* 4 1 

V * 5 J 


M 7 = -cos 


r x. 


K X - , 


M; = a: 45 cos^ 4 


M g 7 = -AT 45 sin^ 4 


M„= - cos 


^* 4 

— * 7 s 

^ 4 


cos^ 4 
sin t ] 4 
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M 9 8 = sin tj 5 + H 5 6 cos rj 5 
Mg = cosrj 5 -H 56 sin t ] 5 

Mg 10 = - sin — 1) 6 

{ X « J 

"i 

Mg 11 = -cos —J) 6 

x* 

\ 6 J 

M 10 — ^ 5,6 COST } 5 
M 10 =~ K 5,6 sin/ 7 5 

A 

M/o 0 =-cos ?J 6 

X* 

V 6 

(x ) 

M\l = sin — ^ 

x* 

V 6 J 

MU =costi 6 
and 

M“ = - sin r) 6 
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The eigenvalues j3 m are the solutions of the transcendental equation obtained by setting 
the determinant of the matrix in Equation 2-47 equal to zero. To obtain the remaining 
coefficients C m of Equation 2-28, we exploit orthogonality and the initial condition. The 

eigenfunctions X jm (x) satisfy the orthogonal relation [Ozisik, 1993] 


'L— i X,Jx)X, Jx)dx 

S — 1 OC ; 

1 — 1 ! Xj 


I" N m ,m = n 
[ 0, m ^ n 


The normalization integral N m and is given by 


(2-48) 


N 


m 



(2-49) 


Wlien t = 0 


0) = ’Lc„x,jx) = -ejx) 


(2-50) 


m = 1 


Now multiplying both sides of Equation 2-50 by —X jn (x), integrating from x t to x, +l , 

a i 

and summing the result from / = 1 to 6 (the number of layers) yields 


6 u x ‘* i 6 b x i+i m =° = 

S-T fy U)X,A*)dx = X— J S C.X i Jx)XJx)dx 

i=i OCj x m= i 


-1 a , , 


6 I , ->.±1 


I C -I— )• ■ (2-51) 

m = 1 ;=1 fy 


: 0 + 0 + ... + C „,N + 0 + 0 + ... 
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Thus 


1 


(2-52) 


C = — 


N„ 


6 K -i+1 

X— \dJx)X m Xx)dx . 

13 a, J 


At tliis point all of the parameters needed to compute 9 r j (x,t) from Equation 2-28 are 
known. The complete solution of the transient problem is given by Equation 2-21. 

The C++ program homog.cpp listed in Appendix A is implemented using the preceding 
results to determine the temperature profile through the thermal radiation detector. A 
parametric study is performed to determine the optimum value of the thickness of the 
actual thermal impedance layer (Parylene). The influence of the contact conductance on 
the accuracy of the results obtained assuming perfect contact between layers is also 
considered. 


2.3 Results 

The actual values of the contact conductances are currently unknown for the fabrication 
technique (sputtering) anticipated for the notional detector. Further studies of the type 
being pursued by Sorensen will be necessary in order to determine these values. The 
models developed in the current chapter are anticipated to be an essential tool of 
Sorensen’s research. The current objectives are to simulate the thermal behavior of the 
detector using different assumed values of the contact conductances to investigate the 
sensitivity of the overall behavior to this parameter. The bulk properties of the materials 
are used in this preliminary study. Assuming continuum mechanics, the thermal 
properties of thin-film layers are known to be a function of the layer thickness [Nath and 
Chopra, 1974; Kelemen, 1976; and Kumar et al, 1994], However, thin-film properties for 
the materials in this investigation are currently unavailable. Tliis is the topic of Chapter 3. 

Figures 4 and 5 show, respectively, the transient temperature distribution of the detector 
with perfect contact and with imperfect contact conductance (in which case the contact 
conductance is arbitrarily set at 50,000 W/m 2 K for the six interfaces). This assumed 
value of contact conductance produces an increase of the temperature at all the locations 
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of the detector, as shown in Figure 5. The temperature is uniform across each junction 
and discontinuous at the interfaces for the non-perfect thermal contact between layers. 

Figure 6 presents the steady-state response of the detector for different values of contact 
resistance as well as for perfect contact. As expected, the lower contact conductance 
results in a higher overall temperature rise at the active junction. It is apparent that the 
contact resistance effect, if present, cannot be neglected without underestimating the 
response of the detector. Also, the curious situation in Figures 5 and 6 in which the two 
components of the thermocouple junction are at different temperatures raises a question 
about how to interpret the temperatures in terms of the thermoelectric effect. 



Figure 4: Temperature response of the notional detector of Figure 1 with perfect thermal 
contact between all layers. 
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Figure 5: Temperature response of the notional detector of Figure 1 with a uniform 
thermal contact resistance between the layers. 



Figure 6: Steady- state temperature response of the notional detector of Figure 1 for 
different values of interlayer contact conductance. 
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The transient temperature response of the active junction, when the detector is subjected 
to a uniform heat flux of 1 W/m on the upper surface of the absorber layer, is shown in 
Figure 7. The response is quasi-first order. The time constant is defined as the time for 
the temperature to undergo 63 percent of the temperature change from the initial to the 
steady-state value. The responsivity of the detector is the difference between the steady- 
state temperature of the active junction and that of the reference junction provoked by the 
heat flux. Figure 8 shows a linear increasing dependence of the responsivity on the 
thickness of the thermal isolation layer. However, the time constant also increases with 
thickness. 



Figure 7: Transient temperature response of the active junction of the notional thermal 
radiation of Figure 1 detector with no contact resistance (thermophysical properties given 
in Table 1). 


The use of less dense materials such as Lare-Si and aerogels, having superior thermal 
insulation properties, is anticipated to be the subject of future investigations. These 
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materials could possibly lead to a lower thermal capacitance for a layer of equal 
thickness, resulting in both a faster time response and a higher responsivity. One of the 
problems with these new materials is that their thin- film thermal properties are unknown. 

The numerical model developed in the current chapter is modified in Chapter 4 in order 
to take into account a nonhomogeneous thermal impedance layer consisting of a mixture 
of aerogel and carbon, or a mixture of Larc-Si and carbon, or a mixture of Parylene and 
carbon. 




■4— ' 
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Figure 8: Time constant and responsivity of the notional detector of Figure 1 and Table 1 
as a function of the thickness of the thermal impedance layer. 
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Chapter 3: The Thin-Film Effect 

3.1 Review of thin-film heat transfer 

In Chapter 2 the effect of contact resistance on detector performance was considered. In 
the current chapter we take up the question of the possible departure of thermophysical 
properties from their bulk values due to the so-called thin-film effect. 

The heat conduction mechanism occurs in solids through collisions among phonons and 
electrons. In metals the main energy carriers are free electrons, while in insulators and 
semi-conductors, the energy is mainly carried by phonons. The macroscale heat transfer 
theories such as heat diffusion given by Fourier’s law fail to describe the heat transfer 
process in micro structures if they are not used appropriately. The classical Fourier’s law 
may be written 


q = -kVT , (3-1) 

where q is the heat flux, k is the thennal conductivity and, VT is the local temperature 
gradient. 

In situations where the characteristic length or/and characteristic time of a structure 
approaches, respectively, the mean free path or the mean free time of the main energy 
carrier, microscale heat transfer models are more appropriate. The classical macroscopic 
Fourier’s law defined by Equation 3-1 breaks down because the definition of thermal 
conductivity depends on the existence of a gradient of temperature within the structure 
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[Majumdar, 1993]. The macroscopic models assume that the physical structure is so large 
and the time scale sufficiently long that sufficient collisions occur between the energy 
carriers as to assure that a local thermodynamic equilibrium is reached. In fact, the very 
concept of temperature at a given point is strictly defined only under the condition of 
local thermodynamic equilibrium. Therefore a meaningful temperature can be defined 
only at points separated by at least the mean free path of the energy carriers [Bejan, 
1988]. We conclude that for heat conduction across a thin-film layer whose thickness 
approaches the order of magnitude of the mean free path, the thermal conductivity of the 
material cannot be defined as given by Equation 3-1. Because an insufficient number of 
heat carriers are present within the material, the temperature field becomes discontinuous 
through the layer and the temperature gradient concept looses its physical meaning within 
the layer. Microscale heat transfer can occur either in length scale or in tune scale or in 
both scales. The two-step phonon-electron interaction model, the phonon- scattering 
model, the phonon radiative model, and the thermal wave model are some of the 
available microscale heat transfer models and are discussed by Tzou [1997], 

An important issue of microscale heat transfer is the dependence of the thermophysical 
properties of thin- films on the micro structure of those films. In a 1984 effort, Decker et 
al. report that the measurements of the thermal conductivities of SiC>2 and AI2O3 thin 
film s are apparently one or two orders of magnitude lower than those for the 
corresponding bulk materials. They also report that the thickness dependence is more 
pronounced for thermal conductivity than for the heat capacity and density in dielectric 
thin- films. Nath et al. [1974] report similar behavior by measuring the thermal 

o 

conductivity of copper films ranging in thickness from 400 to 8000 A m the temperature 
range 100-500 K. This decrease is attributed to a structural disorder inside the thin film 
layer, to a large interface thermal resistance, or to the limitation of the mean free path. 
Oram et al. [1998] also report a decrease of the thermal conductivity of Z 1 -O 2 films with 
thickness ranging from 750 to 10,000 A deposited on AI2O3 substrates. 

Based on kinetic theory, the thermal conductivity of dielectrics and semiconductors is 
given by [Zima, 1960] 
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k = —Cvl 
3 


(3-2) 


where C is the phonon specific heat, v is the speed of sound, and l is the phonon mean 
free path. 

Equation 3-2 shows that the thermal conductivity is directly proportional to the mean free 
path of the mam energy carrier. In thin films the mean free path is reduced from its bulk 
value because of the scattering on the boundaries, and thus so is the thermal conductivity. 
Also the mean free path being smaller in the transverse dimension of the film, the thermal 
conductivity in that direction is smaller than that in the longitudinal direction. For a film 
thickness much larger than the mean free path, the scattering being small, the decrease of 
the thermal conductivity by scattering can be neglected. From Equation 3-2, it is clear 
that the computation of an effective thermal conductivity of the thin film requires the 
knowledge of the mean free path of the main heat carrier. However the contradiction of 
using Equation 3-2 to explain departures from bulk thermal behavior in a mean free path 
regime where Fourier’s law itself probably does not apply is obvious and must be 
recognized. 

Kumar et al. [1994], using the Boltzmann transport theory, derived closed-form 
expressions that predict the reduction in the longitudinal thermal conductivity (in the 
direction parallel to the plane bounding the thin film) of thin metallic films due to 
boundary scattering. Flik and Tien [1990] argue that since the thermal conductivity is 
directly proportional to the mean free path, the reduction of the thermal conductivity 
should be equal to the reduction of the mean free path. They use geometric assumptions 
to evaluate the reduction in the component of the mean free path along the longitudinal 
direction of the film due to the termination of the path lines at the boundaries. They 
derive for the longitudinal conductivity of a thin metallic film, assuming that the path 
lines originate uniformly along the transverse direction (direction normal to the plane 
bounding the thin film) of the thin film. 


]z 

* film 

]z 

* bulk 


1 + — 8 In 
n 


1 + C) + S 
1 + ^-5 


-—cos 1 S-— — (l-s 3 ), 8 < 1, 
n Zk 8 


(3-3) 
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and 


~^- = 1-— - ,S> 1, (3-4) 

hum lx 8 

where 8 is the ratio of the film thickness to the mean free path, and s = (l - S 2 in 
Equation 3-3. 

Figure 9 shows the ratio of the “thin-film” longitudinal conductivity to the bulk 
conductivity for a range of thicknesses determined using Equations 3-3 and 3-4. The 
figure clearly shows the dependence of the thin-film thermal conductivity on the ratio of 
the thin- film thickness to the mean free path length. 



Figure 9: The ratio of the longitudinal thin-film to the bulk conductivities as defined by 
Equations 3-3 and 3-4. 

3.2 The thin-film thermal conductivity anomaly as a roughness effect 

The thermal conductivity is a parameter that has a direct influence on the accuracy of 
models of heat conduction processes. The experimental studies reported in the open 
literature that attempt to recover the effective thermal conductivity of a thin film fail to 
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take into account the roughness of the thin film. They seem to assume that the thin-film is 
completely smooth or that its roughness is negligible. Closed-form equations of the 
thermal conductivity of thin films, such as Equations 3-3 and 3-4, require the knowledge 
of the mean free path of the energy carrier, a property that is not available in the literature 
for the materials of the notional thermal radiation detector. When the film layer is very 
thin (on the order of microns), its relative roughness can no longer be neglected. It may 
be hypothesized that roughness plays an important role in the apparent deviation of the 
properties of individual layers from then bulk values. In other words, reported deviations 
of thermal properties from the bulk values may be an artifact of experimental techniques 
that ignore the possibility of roughness effects. 

The objective here is to determine whether or not unsampled surface roughness could 
possibly explain the observed decrease of thermal conductivity. In fact, no effective 
conductivity data are available in the open literature for the special materials, such as 
Larc-Si and Zinc-antimonide, that are being considered for use in the notional detector. 
One can then imagine performing an experiment to define the statistical description of 
surface roughness for a given sample and then using this information and the bulk 
thermal conductivity with an appropriate model to predict the effective thermal 
conductivity. 

A statistical technique is required for solving the heat transfer equations describing the 
roughness effect for very complicated surface topographies. 

3.2.1 The random walk approach 

The Monte-Carlo method, also called the method of statistical trials, is a system of 
techniques that enables complex physical models of problems to be solved in relatively 
simple manners. The computer is instructed to do most of the work. This approach is 
limited only by the availability of adequate computer resources (mainly speed). In 
Monte-Carlo methods, instead of directly solving an analytical problem, one “plays a 
game” following rules similar to those governing the actual physical process. The game 
has the same outcome as the actual physical process but is in some sense easier to play. 
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The random walk is a Monte- Carlo method used to solve boundary-value and initial 
condition problems. The random walk method enables problems with irregular 
geometries and multiple dimensions to be solved using algorithms that are quite simple 
compared to other methods such as the finite difference or finite element methods. 


3.2.1.1 Presentation of the technique 

In the current study, a two-dimensional steady-state heat conduction problem with no 
internal heat generation is considered. 

The mathematical formulation of the problem is 


d 2 T(x,y ) d 2 T(x,y ) 
a* 2 dy 2 


(3-5) 


with either 


T = f(x,y) (3-6) 

or q(x, y) = 0 (3-7) 

specified on the boundary f. That is, the kind of boundary condition considered is either 
specified boundary temperature or an insulated boundary. 



• internal nodes 
O boundary nodes 


Figure 10: Network of nodes for an arbitrary irregular geometry. 
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One would like to compute the temperature at a given point P(x 0 ,y 0 ) in the two- 

dimensional domain limited by the boundary T, as illustrated in Figure 10. The random 
walk technique works as follow. 

Starting at P(x 0 ,y 0 ) a random number uniformly distributed between zero and unity is 

drawn to determine in which direction to “step”. In a two-dimensional net, one can step 
with equal probability either to the north, the south, the east or the west with a step size 
Ax in the x direction and Ay in the y direction. Having arrived at the appropriate 
neighboring node, the same process is repeated to determine the direction of the next 
step. This process is repeated until a boundary node is reached. Having arrived at a node 
whose location is approximately on the boundary, the next step depends on the type of 
boundary condition: either a specified temperature or an insulated boundary is 
considered. When an insulated boundary node is reached, the next node in the random 
walk is identical to the previous node, as in specular reflection. When a boundary node 
with a specified temperature is reached, a counter ^ 0 is incremented by the value of the 
temperature specified at that node. Then the whole process recommences at point 
P(x 0 ,y 0 ) and is repeated until a large number N of random walks has been completed. 

Finally, after a sufficiently large number of random walks initiated a point P(x (> ,y (; ), the 
temperature at point P(x 0 , y 0 ) is estimated as /N . 

The precision of the temperature computed this way depends on the value of the number 
of random walks N and the fineness of the grid. The larger N , the more accurate the 
result for a sufficiently fine grid. On the other hand increasing either N or the fineness of 
the grid increases the computational time. A trade-off clearly exists between the precision 
sought and the computer tune needed to get that precision, assuming the grid is 
sufficiently fine to provide the deshed precision. 

Another problem of using random walk methods is the lack of reproducibility of the 
results if A is not sufficiently high. When running the same problem several times, one 
never gets exactly the same results. This is due to the difference in the sequences of 
random numbers generated for each experiment. 
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The interested reader is referred to the book by Shreider [1966] for additional information 
on this technique. 

The flowchart given in Figure 11 illustrates the random walk technique in a two- 
dimensional heat condition problem without internal heat generation. 



Figure 11: Flowchart of a two-dimensional random walk (continued on next page). 
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Figure 11 (continued): Flowchart of a two-dimensional random walk. 
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3.2. 1.2 Results 


The program randwalk2d.cpp listed in the Appendix B can be used to determine the 
temperature profile for different geometries. The random walk method developed is 
benchmarked using a problem for which exact closed-form solutions are available. 

The benchmark problem considered using this technique is one-dimensional steady- state 
heat conduction in a flat plate, as illustrated in Figure 12. 



Figure 12: One-dimensional, steady-state heat transfer in a flat plate. 


The analytical temperature distribution is 


WrO^lOOy + 50 


(3-8) 


It is convenient to define an error e e 


1 y \r n Xy i )-T' JcaC '(y i ) 
n y h T emct {y,) 


x 100% 


(3-9) 
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In Equation 3-9, T rw (y) is the value of the temperature obtained using the random walk 
method and n y is the number of discrete divisions in the y direction. 

The random walk code is executed for different values of the number of random walks, 
N . The reproducibility at a given N is studied by running the code several times with 
different starting seeds in the pseudorandom number generator used to draw the required 
random number. The error e is computed for each run. 

Table 2 displays the results obtained for fives runs at a given value of N and shows that 
the error is different from one run to the next. The lack of reproducibility from one run to 
the next is illustrated in Figure 13. 


Table 2: Error values for different runs of the random walk solution to the problem 
defined by Figure 1 1 with n y = 50. 


Number of 
random walks, 
N 

error, 8 (%, Equation 3-9), n y = 50 


first run 

second run 

third run 

fourth run 

fifth run 

50 

19.26 

17.91 

16.91 

16.16 

10.57 

100 

12.48 

7.07 

10.22 

11.61 

10.53 

250 

8.72 

5.71 

8.67 

6.34 

7.03 

500 

5.81 

4.19 

3.33 

4.29 

5.20 

1000 

4.20 

3.78 

2.90 

3.27 

2.76 
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Position y(m) 


Figure 13: First two runs for N = 50 compared with exact temperature distribution for 
the problem defined by Figure 12. 



Position y(m) 


Figure 14: First two runs for N = 1000 compared with the exact solution for the problem 
defined by Figure 12. 
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Table 2 and Figure 14 show that the lack of reproducibility decreases significantly when 
the number of random walks increases. The value of n y is fifty for each of these 
experiments. 

Table 3 presents the influence of the number of random walks on the machine tune (SGI; 
Model: Indigo2 XZ, Extreme; Operating System: IRIX Release 6.5) and on the mean 
error. The machine time is the time it takes to complete the computation, and the mean 
error <£> is obtained by setting T rw in Equation 3-9 equal to the arithmetic mean value of 

T rw , <T, h>, for all the inns; that is, 

1 <T n jy,)> -T exa Ay, 

r exact (y i') 

and < T r Jy i )> = —^T l Jy,) , (3-11) 

n r ,= i 

where n, is the number of runs for a given number of random walks, N. 

As one would anticipate, by increasing the number of random walks, the machine time 
increases and the result becomes more accurate. Once again the number of y divisions, 
%, is fifty for each of these numerical experiments. 


<£> 


i=l 


x 100% 


(3-10) 


Table 3: Influence of the number of random walks on the machine time and the mean 
error for the problem defined by Figure 12 with n y = 50. 


Number of random 
walks, N 

Machine time (s) 

Number of runs, n, 

Mean error (%) 

50 

8 

5 

7.88 

100 

17 

5 

3.61 

250 

44 

5 

2.90 

500 

88 

5 

2.52 

1000 

176 

5 

1.73 

5000 

885 

1 

1.54 


The Thin-Film Effect 


41 

































Figure 15 shows a comparison between the result obtained for just one ran with 5000 
random walks and the exact solution. The error being only about 1.5 percent, the 
temperature obtained by the random walk technique can be acceptably accurate. 
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Figure 15: Comparison of the temperature distribution computed using the random walk 
method with N = 5000 and n y = 50 with the exact solution for the problem defined by 
Figure 11 with. 
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The two-dimensional random walk technique has also been nsed to compute the 
temperature distribution when surfaces are modeled as trigonometric functions, as 
illustrated in Figure 16. 

Figure 17 represents the temperature profile at x = 5 pm in Figure 16 computed using the 
random walk method for different values of the number of walks N. 
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Figure 16: Irregular geometry representing a layer of variable thickness. 



Figure 17: Comparison of the temperature distribution computed using the random walk 
method at x = 5 pm in Figure 16 with n y =50, N = 1000, 5000, and 10000. 
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3.2.2 Effective thermal conductivity of a thin-film layer of irregular thickness. 

The investigation of the effective thermal conductivity of a thin-film layer based on a 
roughness effect begins with the simple two-dimensional model illustrated in Figure 18. 
The idea is that if the results obtained using this model are consistent with the expected 
trend, that is a decrease of the thermal conductivity of the thin-film as the roughness 
increases, then more realistic surfaces would be considered. It is convenient to use a finite 
difference model in order to determine the temperature distribution for the geometry in 
Figure 18. The finite difference method is conceptually simple for problems having 
regular rectangular boundaries. It involves the use of nodal networks, finite difference 
approximations for derivatives in space and tune, standard energy conservation 
formulation concepts, and computer solution of systems of algebraic nodal equations. 
The control volume approach is used for the discretization of the governing equations. 
The C++ code ejf cond.cpp listed in Appendix C and developed for computing the 
temperature profile using the finite difference method in a two-dimensional regime 
without internal heat generation is adapted from Vick [19981. If this study indicates that 
an unsampled roughness effect can explain the observed thin-film effect, then the random 
walk method will be used to study more realistic geometries. 


T=100 K 
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Figure 18: Model of a hypothetical two-dimensional rough surface. 
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Here, the effective conductivity k eff of a real surface is defined as the conductivity of a 
hypothetical plane surface having the following three properties: 

• the same length L as the surface in Figure 18, 

• the same total heat flux as the problem defined by Figure 1 8, 

• a uniform thickness of H + R / 2, where H and R are defined in Figure 18. 

A discretized version of Figure 18 is shown in Figure 19. Because the problem is two- 
dimensional, the total heat flux in the y direction is appro ximated using the bottom row of 
control volumes in the Figure 19. Then for a given control volume at the bottom of 
Figure 19, the heat flux through the ith heat flow channel is appro ximated by 


Qi 


T(Xj,Ay 12)-^ 


v bulk 


Ay/2 


-Ax 


(W/m). 


(3-12) 


The total heat flux is then approximated by 

Qt = x<2, = 2^f j [{T(x i ,Ay/2)-T l )Ax,] (W/m), (3-13) 

i Ay r 


where ii, is the number of control volumes in the x direction, and k bM is the bulk 
conductivity of the material. 

The total heat flux of the hypothetical flat surface is given by 


Q — k eff L 


T -T 
1 2 ' 1 

(H +R/2) 


(W/m) . 


(3-14) 


Equating Equations 3-13 and 3-14 yields 



(H + R/2) 2 

L T 2 -T y 


£[(r(.r,,Ay/2)— TjAr,] 


(3-15) 
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Insulated 



Figure 19: Control volume discretization of the hypothetical two-dimensional rough 
surface of Figure 18. 



R/(R+H) 


Figure 20: ratio of the effective thermal conductivity to the bulk thermal conductivity for 
the problem defined by Figures 18 and 19 with forty control volumes in the x direction 
and LIH = 5. 
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Ill order to study how the roughness, represented by the nondimensional quantity 
R/(R + H) influences the ratio k eff /k bM , the value of R is set constant in the finite 

difference procedure and the ratio is computed for different values of H . Figure 20 
illustrates the result obtained. 

The trend of Figure 20 is completely the opposite of what is expected. That is if 
unsampled roughness is to account for the observed [Nath and Chopra, 1974; and Orain 
et al., 1998] behavior of effective thermal conductivity in thin- film layers, a decrease of 
the effective thermal conductivity is expected as the roughness increases. The effective 
conductivity for the model of Figure 19 exceeds the bulk value when the roughness 
increases, and approaches the bulk value as the surface becomes smoother . Looking 
more closely at the geometry Figure 18, one can see that when H becomes smaller 
compared to R, the surface at higher temperature 77 is closer to the plane at temperature 
T 2 so that the local gradient of temperature is larger. The larger temperature gradient 
yields a higher total heat flux and therefore a higher effective conductivity. 

We conclude that the approach used is unable to reproduce the trends reported in the 
literature. In other words, unsampled surface roughness cannot explain experimental 
results. Therefore, the random walk technique is not considered further for studies with 
more realistic surface models. 

Evidently other approaches such as experimental measurements or parameter estimation 
techniques will have to be used to determine how the thermal conductivities for the 
materials of current interest are affected by the thin- film effect. 

In Chapter 4 we investigate the effect on effective thermal and electrical conductivity of 
doping the thermal impedance layer with conducting material such as finely divided 
graphite. 
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Chapter 4: Effect of Doping the Thermal Impedance 

Layer with Conducting Particles 


In the previous Chapter we investigated and rejected surface roughness as a possible 
explanation for the so-called “thin-fihn” effect. In this Chapter we consider the effect on 
thermal and electrical conductivity of doping the thermal impedance layer with finely 
divided carbon. 

The NASA team led by Mr. Edward H. Kist, Jr., has proposed that the thermal impedance 
layer of the thermal radiation detector be doped with up to fifty percent by volume of 
finely divided graphite. The reason for this is to decrease the internal electrical resistance 
of the thermal impedance layer, thereby assuring the electrical connection between the 
two junctions of the thermal radiation detector. As stated before, the geometry and size of 
the detector preclude linking the two junctions by a more traditional conductor. 
Experiments carried out by Sorensen [1998] showed a decrease of the electrical 
resistance of thin layers of Larc-Si doped with powdered graphite. Minimization of the 
electrical resistance of this layer is a key factor for reducing an important noise source in 
radiation detectors. Johnson noise in an electric circuit is known to be directly 
proportional to the square root of its electrical resistance [Lenoble, 1993]. The noise 
sources in thermal radiation detectors are discussed by Lenoble [1993] and Haeffelin 
[1997], 

An obvious concern is that reduction of the electrical resistance would also lead to 
reduction of the thermal resistance and concomitant reduction of the detector’s 
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responsivity. The objective here then is to develop an analytical model that predicts the 
performance of the detector whose thermal impedance layer is doped with conductive 
particles. This situation is shown in Figure 21. The model is applied to the specific case 
of a Parylene impedance layer doped with carbon particles. 



■ Absorber layer (10 pm) 

■ Zinc-antimonide layer (1 pm) 

■ Platinum layer (1 pm) 

m Doped thermal impedance layer 
(25 pm) 


Figure 21: Detector geometry with graphite doped thermal impedance layer. 


The behavior of the thermal radiation detector is studied first using a model of the 
detector derived from a model by Vick and Scott [1998] and then using micromechanical 
models. 


4.1 Heat transfer in a heterogeneous material modeled as a coupled, two-step 
process. 

4.1.1 Formulation 

The transient thermal response of a heterogeneous material depends on the micro structure 
of the material. The one-dimensional heat transfer model developed by Vick and Scott 
considers relatively small spherical particles (the Biot number of less than 0.1) embedded 
in a matrix. In this model heat is first conducted to the matrix material from the 
boundaries and then transferred from the matrix to the particles through a contact 
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conductance. The particles are not in perfect thermal contact with the matrix material and 
they can receive heat only from the matrix. The boundaries of the thermal impedance 
layer in the current version of the model are the two zinc-antimonide layers of the active 
and reference junctions shown in Figure 21. 

On the basis of an energy balance Vick and Scott have proposed that 

C„^ L = K^ L -H(T m -T„) (4-1) 

at dx 

and 

C P = H (T m — T p ) (4-2) 

at 

for the doped matrix for one-dimensional transient heat flow. 

In Equations 4-1 and 4-2 C m = (1 - f)(pc) m is the heat capacity of the matrix per total 
volume and K = (1 - f)k m is the thermal conductivity of the matrix-particle system. 
C p = f(pc) p is the heat capacity of the embedded particles per total volume, and /is the 

volume fraction of the embedded particles. The coupling coefficient between the matrix 
and the particles is H = h(NA p ) , where A p = n a 2 is the surface area of a particle, h is the 

contact conductance, and N is the number of particles per volume. 

This model is analogous to the two-step model developed by Tzou [1997] in order to 
describe the micro scale heat transfer process for metals in which the electron gas is 
heated first and then the heat is transferred to the metal lattice. It holds for each of the six 
layers of the thermal radiation detector, but the number of particles is set equal to zero for 
layers with no carbon particles. The boundary conditions are a specified temperature at 
the base of the detector and a constant and uniform heat flux applied at the upper surface. 

The coupled equations, Equations 4-1 and 4-2 are not solvable analytically. Vick and 
Scott used a fully implicit finite difference numerical approach to discretize those 
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equations. The fully implicit discretization makes the results unconditionally stable. The 
results of the discretization are 


(T . 
C_ Ax- m ’ J 


-T’~ l ) 

m, / / 


At 


/C 


r' 

J mj+l 


-T } 

m,j 


Ax 


T’ 

-K- m ’ ] 


-T' 

1 m,j - 1 


Ax 


-^(7’;,, -r;,.) 


(4-3) 


and 


-i 


C P Ax 


T h ~ T L 

At 


= HAx{T' mJ -T' PJ ) 


(4-4) 


where j and i are the indices of the nodes and time steps, respectively. 

To solve these coupled equations, T Pj from Equation 4-4 is first substituted into 
Equation 4-3 and the resulting tridiagonal matrix is solved for the temperatures T’ mJ at 
each time step. The values of T’ n j are then substituted into Equation 4-4 to compute the 
particle temperature T P ] within the thermal impedance layer. 

The program heterog.cpp listed in Appendix D determines the temperature response of 
the detector for different values of the contact conductance between the carbon particles 
and the matrix material. The values of the contact conductance and the particle diameters 
are chosen so that the Biot number always remains small (less than 0.1). In this case the 
lumped capacity approximation is valid within each particle; that is, temperature 
gradients are negligible within a given particle. 

4.1.2 Results 

In all the cases considered in this chapter, the contact resistances between the different 
layers of the thermal radiation detector are assumed negligible, and the carbon particles 
are assumed to be uniform in size. It is recognized that this latter assumption would be 
difficult to realize in practice. However, the assumption approximates the situation in 
which a mean particle size is used to represent some size distribution. The contact 
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resistances between the layers forming the detector are neglected in order to focus on 
only one phenomenon at a time. 


Studies of the effect of three parameters on detector response were carried out using the 
two-step model process. The three parameters are: 

• the volume fraction, f, of the particles, 

• the contact conductance between the particles and the matrix and 

• the particle size. 



2 

Figure 22: Temperature response of the active junction with h = 50 W/m K. 


The model developed here holds only for a matrix with dimensions much greater than a 
particle diameter and with a small volume fraction of particles within the matrix. Vick 
and Scott do not specify the volume fraction limit at which the model breaks down. 
Therefore the model is used here only for volume fractions of particles ranging between 
zero and thirty percent. 
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Figure 22 shows the transient temperature response of the active junction for four cases 
of particle volume fraction, f. The contact conductance between the particles and the 
matrix is 50 W/nr K in all four cases. The figure reveals an increase of the steady- state 
temperature of the active junction of the detector as the particle volume fraction 
increases. This result seems to be counterintuitive; that is, it is expected that the thermal 
impedance would decrease as the volume fraction of conductive particles increased. 



Figure 23: Effect of contact conductance between the particles and the matrix on the 
steady- state temperature rise of the active junction. 


Figure 23 shows that the contact conductance between the matrix and the particles has 
negligible effect on the response of the thermal radiation detector. The greatest change of 
the response of the detector with respect to the contact conductance is less than one 
percent for the range of particle volume fractions considered in this study. 
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For the geometry and range of parameters investigated, the contact conductance has more 
effect on the temperature of the particles than on the temperature of the active junction of 
the detector, as shown in Figure 24. The heat is transferred from the matrix to the 
particles; therefore when the contact conductance increases, the time lag of the particles 
decreases. Note that the curve labeled “matrix” in Figure 24 actually consists of three 
superposed curves corresponding to the three values of h. Therefore it is clear that the 
contact conductance has little effect on the matrix temperature. The time lag of the 
particles decreases rapidly as the contact conductance increases. When h is about 10 W / 
m K the tune lag between the matrix and particle temperatures is effectively zero. 



Figure 24: Matrix and particle temperatures at x = 12 pm for different values of the 
contact conductance. 


The effect of the particle size on the steady-state response of the detector is shown in 
Figure 25. The range of particle diameters was chosen so that the lumped capacity 
assumption holds. In this range it was found that the particle size has a very small effect 
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on the output of the detector. The maximum temperature change of the active junction is 
less than one-third of a percent. 



Figure 25: Particle size effect on the steady- state temperature rise of the active junction 
for a volume fraction of particles of ten percent. 


4.1.3 Discussion of the results 

The temperature response shown in the Figure 22 for the thermal radiation detector does 
not agree with what one would expect by doping a material of low thermal conductivity 
(the conductivity of the Parylene is 0.084 W/m K) with a material of higher thermal 
conductivity (the conductivity of the carbon particles is 1.59 W/m 2 K). In Equation 4-1 of 
the two-step transient model, the conductivity K is equal to (1 - f)k m , where k m is the 

conductivity of the matrix. This expression for K requires that the thermal conductivity 
decrease as the volume fraction of particles increases, with the result that the temperature 
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of the active junction increases. With this model, the steady-state temperature distribution 
depends only on the thermal conductivity of the matrix and the volume fraction of 
particles, and not on the thermal conductivity of those particles. The steady-state 
temperature increases as the volume fraction of the particles increases. This is contrary to 
intuitive expectations as well as to observations of decreases in electrical resistance of 
Larc-Si layer with an increased loading of graphite particles [Sorensen, 1998]. 

It is hypothesized that this two-step model could be unproved by somehow taking into 
account the thermal conductivity of the particles, even if the lumped capacity 
approximation is assumed. One way would be to solve the Laplace equation for heat 
conduction subject to the conditions of continuity of the temperature (assuming perfect 
contact between the particles and the matrix) and heat flnx at the interface of each 
particle and the surrounding matrix. This direct method would of course be 
computationally intensive. The alternative approach in this thesis, described below, is to 
compute equivalent properties for a heterogeneous thermal impedance layer. The 
response of the detector is then analyzed using the analytical transient model developed 
in Chapter 2. 


4.2 Heat transfer in a heterogeneous material using micromechanical property 
models. 

4.2.1 Thermal analysis 

The determination of the thermophysical properties of heterogeneous mixtures has been 
and continues to be a very active area of research because of the wide range of 
applications of these materials in engineering [Meredith et al., 1962], Although some 
properties of these mixtures, such as the specific heat, do not present any difficulties, 
other thermal and transport properties, such as thermal or electrical conductivities, 
require more sophisticated treatments. The simple mixing rales that consist of averaging 
the transport properties of pure phases to get those of the mixture do not apply in these 
latter cases. The complexity of the problem increases for more concentrated mixtures. 
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The primary interest here is the determination of the effective thermal conductivity and 
heat capacity of a two-phase composite material consisting of spherical carbon particles 
embedded in a continuous matrix. 


4.2.1.1 Effective thermal capacity of the mixture 

The effective heat capacity of a mixture can be directly obtained using a simple mixing 
rule. The effective specific heat of the matrix-particle mixture is 

~c = f' P c P + f' m c m (4-5) 

where subscripts p and m indicate “particle” and “matrix” respectively and f' is the 
appropriate mass fraction, 



IYI 

mixture 


Alternatively, we can define the volume fraction 
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V, 


v„ 


(4-6) 


(4-7) 


with i = p or m. 


The relation between f' and/) is f'= -=^f, 

P 


where p is the density of the mixture. 


- p P v P +p m v m 

P = ~ L — £ 

v 

mixture 


(4-8) 


(4-9) 
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Equation 4-5 is a statement of conservation of energy, and Equation 4-9 is a statement of 
conservation of mass. 

Combining equations 4-5 and 4-6 yields 

Jc = f{pc) p +(l-f){pc) m . (4-10) 

In Equation 4-10, pc is the effective heat capacity of the matrix-particle mixture and /is 
the volume fraction of the particles. Equation 4-10 is valid for both dilute and 
concentrated mixtures. 


4.2. 1.2 Effective thermal conductivity 

The theoretical investigation of the equivalent thermal conductivity requires the 
distinction between dilute and concentrated mixtures. 


4.2. 1.2.1 Dilute mixtures 

Dilute mixtures are mixtures in which the disturbance of the lines of heat flow by one 
particle is independent of the disturbance provoked by any other particle. For these dilute 
mixtures, Maxwell [1954] derived an equation for the effective electrical conductivity by 
solving the Laplace equation for a single sphere in a matrix. The mathematical behavior 
of the physical properties thermal and electrical conductivity being analogous, the result 
obtained for the effective electrical conductivity is applicable for the effective thermal 
conductivity. 

The result obtained by Maxwell when the thermal conductivity is substituted for 
electrical conductivity, is 


-_ K d +2-2f(\-K d ) 
K d +2 + f(l-K d ) 


(4-11) 
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where K is the relative effective thermal conductivity that is the ratio of the effective 
thermal conductivity of the mixture, k , to that of the matrix , k m . The quantity K<j is the 
ratio of the thermal conductivity of the spherical particles, k d , to that of the surrounding 
matrix. 


Equation 4-11 is valid regardless of the size (particles of uniform size) and spatial 
distribution (random or ordered dispersions) of spherical particles as long as the volume 
fraction of the particles is less than 0.2 and the contact resistance between the particles 
and the matrix is negligible. 

Hasselman and Johnson [1987] modified the theory of Maxwell and derived an 
expression for the effective thermal conductivity of mixtures with a thermal contact 
resistance between the particles and the matrix. They found that the effective thermal 
conductivity of a given composite system and dispersed phase do not depend only on the 
volume fraction of the particles but also on the particle size. They show that 


K d + 
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2k d ^ 
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1 ~K d + 
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K d + 


IK 
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+ 2 


+ f 


1 ~K d + 


ah 


(4-12) 


where h is the contact conductance between the spherical particles and the matrix 
material, k d is the thermal conductivity of the particles, and a is the diameter of the 
particles. Like the Maxwell expression, Equation 4-12 is valid only for dilute mixtures. 
When h goes to infinity Equation 4-12 reduces to Equation 4-11. 


4.2. 1.2.2 Concentrated mixtures 

The determination of the effective properties of composites in winch the volume fraction 
of the dispersed phase is high (f > 0.2) presents increased mathematical difficulties 
[Meredith et al., 1962], As the volume fraction of particles increases, the heat flow fields 


Effect of Doping the Thermal Impedance 
Layer with Conducting Particles 


59 



surrounding the particles begin to interact. The micro structure of the composite material 
begins to play a more important role and so does the manufacturing technique upon 
which the micro structure depends. The effective properties depend on the arrangement of 
the particles with respect to the direction of the heat flow hues. Many predictive formulas 
of the effective thermal conductivity have been derived for different geometries and 
arrangements of the dispersed phase. Lord Raleigh [1964] found that for a cubic array of 
spheres where the heat flow field is on the average perpendicular to a side of the cube, the 
effective thermal conductivity is given by 


K = 1- 


3 / 

(2 + K d )/(l-K d ) + f- 0.525(1 - K d )f 10/3 /(4 / 3 + K d ) 


(4-13) 


Equation 4-13 is valid for / less than n/6. Meredith and Tobias [1962] derived another 
equation for K for the same conditions using a different potential function. They found 
that 


K 1 (2 + K d )/(l - K d ) + f K d )f 1013 /[4/3 + Kj + 0.409(1 - )/ 7/3 ] (414) 

Equations 4-13 and 4-14 are valid for the cases of ordered arrangements of the particles 
where the position of each particle is known and the locations of the particles are chosen 
so that they are symmetric and thus produce fewer interactions among particles. The 
general results for randomly distributed particles are quite different. A large number of 
relations exist in the literature for the computation of the effective thermal conductivity in 
the case of randomly distributed spherical particles in a matrix. The relation derived by 
Meredith and Tobias [1962] is used in the current investigation. For two-phase mixtures, 
Meredith and Tobias found that 


'2(K d +2)+2(K d -\)f~ 

'{2-flK t +2)+2(K l -\)f 

_ 2(K t +2)-(K J -1)/ _ 

_(2-flK d +2)-(K d -\)f _ 


(4-15) 


The maximum and minimum values of the thermal conductivity are obtained in the 
folio whig limiting cases [McLachlan, 1990] : 
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• the particles of the suspended phase link up parallel to the direction of the mean 
heat flow, as in Figure 26 (a). 

• the particles of the suspended phase link up perpendicular to the direction of the 
mean heat flow as in. Figure 26 (b). 




(b) 

Figure 26: (a) lamina parallel to the heat flow direction; (b) lamina perpendicular to the 
heat flow direction. 


The effective thermal conductivities corresponding to these two cases are, respectively, 
the arithmetic mean and the harmonic mean of the thermal conductivity of the matrix and 
the particles. That is. 
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k=fk d +(l-f)k m 


(4-16) 


and 


l/k = f/k d +(l~f)/k m . (4-17) 

In Equations 4-16 and 4-17 k is the effective thermal conductivity, k<j the thennal 
conductivity of the particles and k m the thennal conductivity of the matrix. Note that 
Equations 4-16 and 4-17 follow directly from the limiting cases of the conductance of 
parallel and series networks, as shown if Figure 26. 


4.2. 1.3 Results 

The effective thennal conductivity of dilute mixtures is governed by the magnitude of the 
parameter k d /(ah ) , as described by Equation 4-12. In a matrix doped with very small 
particles this tenn can have a strong influence on the effective thennal conductivity, as 
shown in Figure 27. The effective thennal conductivity is maximum when the contact 
between the matrix and the particles is perfect. The effective thennal conductivity is 
lowest when the nondimensional parameter k rl /(ah) goes to infinity because of the 

additional thennal resistance to the heat transfer that builds around the particles. The 
determination of the contact resistance between the particles and the matrix would be a 
challenging experimental problem because of the dimensions of the detector. To obtain a 
first approximation of the effective thermal conductivity, the contact conductance tenn is 
neglected in this current effort and thus the Equations 4-11 and 4-15 are used in 
determining the effective thennal conductivities. The use of these two equations assumes 
randomly distributed particles in the matrix. The computer code homog.cpp listed in 
Appendix A has been modified to take into account the dependence of the heat capacity 
and the thennal conductivity of the thennal impedance layer on the volume fraction of 
the particles. 
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Figure 27: Effective thermal conductivity of the doped thermal impedance layer for K<j = 
km/kd = 18.93 and for a range of kd/ah. 


The effective thermal conductivity of the thermal impedance layer is shown in Figure 28 
for values of the particle volume fraction of up to fifty percent. An increase of the 
effective thermal conductivity of the matrix occurs as the volume fraction of the particles 
increases. Also presented in that figure are the two limiting cases, Equations 4-16 and 4- 
17. The effective thermal conductivity is bounded by these two limiting cases for the 
entire volume fraction range. It is clear that using one of these two limiting cases would 
lead to important errors in the value of the effective thermal conductivity, assuming the 
other models are correct. Equation 4-15 being in good agreement with the Maxwell 
expression, Equation 4-11, for a volume fraction less than twenty percent, it is used for 
computing the effective thermal conductivity for the entire range of values of the particle 
volume fraction. 
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Figure 28: Prediction of the effective thermal conductivity of the doped impedance layer 
for k m = 0.084 W/m K and k d = 1 .59 W/m K. 



Figure 29: Active junction transient temperature rise of the detector with doped thermal 
impedance using homogeneous approximations method for a range of values of the 
volume fraction of the particles and kd/(ah) = 0. 
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Figure 29 shows that the temperature of the active junction drops as the volume fraction 
of the particles increases, as expected. As the effective thermal conductivity of the 
impedance layer increases, more heat is transferred through the detector. The matrix 
being doped (perfect contact) with better thermally conducting particles, its ability to 
transfer heat to the layers below is unproved. This is in keeping with the electrical 
conductivity experiments reported by Sorensen. Although the heat capacity of the matrix 
increases as the volume fraction of particles increases, Figure 29 shows a faster tune 
response of the detector with increased particle volume fraction; that is the lower steady- 
state temperature is reached more quickly. 

Figure 30 shows how the steady-state temperature of the active junction and the 
corresponding potential difference decrease when the particle volume fraction is 
increased. The temperature drops from about 300 to 75 pK and the potential difference 
from about 0.28 to 0.008 pV as the volume fraction goes from zero to 50 percent. The 
doping of the matrix by the particles greatly decreases the responsivity of the detector. 



Figure 30: Active junction steady-state temperature rise and corresponding potential 
difference (Seebeck coefficient x Temperature Rise) for a range of values of the volume 
fraction of the particles. 
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Filially the temperature distribution through the detector with the thermal impedance 
layer doped with fifty percent of particles is shown in Figure 31. The temperature 
gradients are smaller than those in Figure 4 but the temperatures of the two layers 
forming each junction are still essentially uniform through the junction. While the 
performance of the thermal impedance layer is diminished, the temperature of the 
reference junction remains that of the aluminum substrate. At this concentration of 
particles, the steady-state condition is reached in about ten milliseconds. That is, the tune 
response of the doped detector is almost three tunes faster than that of the configuration 
without particles. 

Figure 32 shows the steady- state temperature distribution through the detector computed 
using the two-step process model and the macromechanical model with thirty percent of 
particle volume fraction, and that of the undoped thermal detector. The contact 
conductance h is fixed high enough so that it can be assumed perfect contact between the 
matrix and the particles. 

The use of the mixing rules for computing the effective thermal capacity and predictive 
formulas for the effective thermal conductivity of the matrix and particles gives results 
that are in accord with intuition. Although the predictive correlations start to break down 
for higher values of volume fraction, results given by the homogeneous approximations 
encompassed by Equations 4-10 and 4-15 are more realistic than those for the two-step 
model reported in Section 4-1. 
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Figure 31: Temperature profile through the detector doped with fifty percent of carbon 
panicles using the homogeneous approximations method to estimate the thermal 
conductivity and the heat capacity of the thermal impedance layer. 



Figure 32: Comparison of steady-state temperature profile through the detector obtained 
using (a) the two-step model, with (c) that obtained using the micromechanical model, 
and (b) the undoped detector model (the thermal impedance layer of the detector is doped 
with thirty percent of particles for the first two cases). 
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4.2.2 Electrical analysis 


4.2.2. 1 Effective properties 

The effective electrical conductivity of the system composed of a matrix doped with 
conductive particles is obtained using the same equations derived for the effective 
thermal conductivity (Equations 4-10 through 4-17). The objective is to determine how 
the effective electrical resistance between the upper surface (Si) and the lower surface 
(S 2 ) of the thermal impedance layer, illustrated in Figure 33, changes as the volume 
fraction of the particles increases. The electric current flow is from Si to S 2 and the other 
two boundaries of the thennal impedance are assumed to be electrically insulated. 


Si 



Figure 33: Thennal impedance layer doped with particles. 


The heterogeneous matrix-particle system is electrically modeled as a homogenous 
material with effective electrical properties. Based upon the effective themal conductivity 
equations, the relative effective electrical conductance is given by 


'2(A„+2)+2(A,-l)/' 

'(2-/XA i+ 2)+2(A„-l)/' 

_2(A,+2)-(A J -l)/_ 

. (2-/)(A rf +2)-(A J -l)/_ 


(4-18) 
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where A is the relative effective electrical conductivity; that is, the ratio of the effective 
electrical conductivity of the mixture, <7, to that of the matrix, o m , and the quantity A d 
is the ratio of the electrical conductivity of the spherical particles, <7 , , to that of the 
surrounding matrix. 

The electrical conductivity for the two limiting cases, i.e. Case (a) and (b) of Figure 26, 
is, respectively 

ct = /o-„+(1-/K, (4-19) 


and 


l/<T = //CT J +(l-/)/<T, 


(4-20) 


The effective electrical resistance (/?) between the surfaces Si and S 2 of the thermal 
impedance layer is a function of the cross-sectional area S, the length L, and the effective 


resistivity 


- 1 
< P = = 
<7 


of that layer; that is, 


R = cp 


w 


(4-21) 


4.2.2.2 Results and Discussions 

The effective electrical resistance is determined using the homogeneous approximations 
for a Parylene layer with embedded carbon particles. The experimental determination of 
the electrical resistance actually underway at NASA Langley Research Center is done for 
a Larc-Si matrix with embedded carbon particles. The thermophysical properties of tins 
new material being currently not well known, Parylene is considered in this investigation. 
However, the same equations could be used latter on when the properties of the Larc-Si 
are known. 
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The resistivities of Parylene and carbon are respectively 6xl0 14 and 7.27xl0 4 Q'm' . 
The resistivity of Parylene being several orders of magnitude higher than that of the 
carbon, Equations 4-18 through 4-20 can be simplified and become respectively, 


<T = AX<7„ 


(1 + /X2 + /) 
(1-/X2-/) 


X ( 7 „ 


(4-22) 


o' = f° d 


(4-23) 


and 


a - 


1 -/ 


(4-24) 


Figures 34 and 35 present the behavior of the effective resistivity and electrical resistance 
of the thermal impedance layer as the volume fraction of the particles increases. Although 
the figures show decreasing resistivity and electrical resistance, these two properties still 
remain far too high. The electrical conductivity of the Parylene is so low (o m = 
0.167 xlO _14 Q _1 m _1 ) [Parylene Coating Services, 1999] that doping it with the carbon is 
not able to sufficiently improve the electrical properties of the thermal impedance. 

Another factor that influences the electrical resistance is the ratio ^ . This multiplying 

factor of the resistivity is equal to 7xl0 3 for the dimensions of the thermal impedance 
(square cross-section with L = 60 jum and w = 25.4 jum). The results show that doping 
the thermal impedance leads to a decrease in the electrical resistance. The decrease would 
be greater if the matrix had a much higher electrical conductivity to begin with. 

The limiting cases, Equations 4-23 and 4-24, are not shown in Figure 31 because of their 
order of magnitude compare to that of Equation 4-22. 

The approach presented here could be used to characterize the electrical performance of 
other combinations of matrices and particles, and different dimensions of the thermal 
impedance layer. 
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Figure 34: Effective resistivity of the thermal impedance layer defined by Figure 33 with 
o m = 0.167 xlO 14 Q -1 m" 1 . 



Figure 35: Effective electrical resistance of the thermal impedance layer defined by 
Figure 33 with L = 25.4 pm and S = 60x60 pm 2 . 
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Chapter 5: Conclusions and Recommendations 


5.1 Conclusions 

Several results follow and conclusions can be drawn from the research reported in this 
thesis: 

1. An analytical model based on the orthogonal expansion technique has been 
developed to characterize the effect on the performance of the thermal 
radiation detector of hypothetical contact resistance between the layers of 
detector. 

2. The model shows that the existence of a contact resistance yields a greater 
temperature response of the detector. For example the steady- state 
temperature response of the detector rises from about 300 pK for a perfect 
contact between layers to 600 pK for an hypothetical 10,000 W/m 2 K contact 
conductance between the layers. 

3. The responsivity and time response of the detector both increase as the 
thickness of the thermal impedance layer increases. This suggests that an 
optimum configuration may exist within the constraints of geometry and input 
signal strength. 

4. In order to investigate the possible effect of interface roughness on thermal 
conductivity, the random walk technique has been adapted for solving 
conduction heat transfer problems in two-dimensional geometries with 
specified boundary temperatures or insulated boundaries. The reproducibility 
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of the method is benchmarked by solving the special case of a one- 
dimensional heat conduction problem in a flat plate. Both the precision of the 
method and the machine time increase as the number of random walks 
increases. 

5. The roughness effect model is unable to reproduce the decrease of the 
effective thermal conductivity of thhi-fOms layers reported in the literature, 
and must be rejected as a possible explanation of the so-called thin- film effect. 

6. A two-step model is used to characterize the heat conduction in the thermal 
impedance layer doped with conducting particles. The model predicts that by 
doping the thermal impedance layer with particles with higher thermal 
conductivity, the steady- state temperature response of the detector would 
increase, which is a counterintuitive. 

7. The detector with doped thermal impedance layer performance is also 
investigated by determining effective electro- and thermophysical properties 
of the matrix-particle system using micromechanical models. Doping the 
thermal impedance layer of the matrix by conducting particles is shown to 
decrease the thermal impedance of the detector, as expected. As a result the 
detector responsivity decreases significantly. However the time response of 
the doped detector is found to be almost three times faster than that of the 
configuration without particles. 

8. The decrease of the electrical resistance of the Parylene layer doped up to 50 
percent by carbon particles does not have a significant effect on the electrical 
performance of the notional detector. 

9. The thermal contact resistance and the micromechanical mixture models are 
expected to be useful in the ongoing parameter estimation activity in the 
Thermal Radiation Group. 
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5.2 Recommendations 


Several recommendations can be drawn from the research reported in this thesis: 

1 For further research in the development of the notional thin-film 
thermal radiation detector that is the subject of this effort, it is 
suggested that a special effort be dedicated to the determination of the 
unknown properties of the materials of the notional detector. 

2 Although the current research was done using macroscale theories, the 
knowledge of properties such as the ratio of the thickness of a layer 
and the mean free path would give more insight in the understanding 
of the heat transfer process in the thermal radiation detector. 

3 Measurement of the contact resistance between layers or the 
integration of its effect in an equivalent thermal resistance of layers 
will also be valuable. 

4 Efforts should also be oriented toward the determination of other 
combinations of materials for the thermal impedance layer in order to 
decrease the electrical resistance of this layer without sacrificing its 
thermal performance. 
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Appendix A: homog.cpp 


The program homog.cpp computes the temperature rise above the aluminum substrate within the 
thermal radiation at a given location, x, and a given time, t. 

List of functions 

* eigen 

The function eigen returns the determinant of the matrix in Equation 2-47. 

* zbrak 

This function is used by the bisection method for root finding. When a function has multiple roots 
in an interval, zbrak determines the subintervals in that interval containing just one root. 

*findroot 

findroot is a root finder by bisection method. This function is used to determine the eigenvalues by 
setting the determinant returned by the function eigen equal to zero. 

* cprime 

cprime computes sum{x_node[j-l ]*(1/K[j-1]-1/K[j])} for j = 1: i. 

* Nm 

Nm is the normalization integral. 

* Cm 

Cm returns the coefficients of 6,j. 

List of variables 

A matrix in Equation 2-47 

A[] and B[] vectors of the coefficients of the X i/n ’s 

alpha [ ] thermal diffusivity vector, m 2 /s 
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beta 

eigenvalue 

guess beta 

guess for the eigenvalue beta 

kij [ ] 

vector defined by Equation 2-41 

K[] 

thermal conductivity vector, W/m.K 

layer 

index of a layer 

Mass density 

mass density vector. Kg/m 3 

nb 

number of eigenvalues 

Q 

heat flux, W/m 2 

Specific heat [ ] 

specific heat vector, J/Kg.K 

t 

time, s 

t trans 

transient component of the temperature rise, K 

t ss 

steady- state component of the temperature rise, K 

Thickness [ ] 

thickness vector, m 

x node [ ] 

layer upper surface coordinate vector, m 

X 

location, m 




* Headers and Library functions 
#define WANT_STREAM 
#define WANT_MATH 

#include "newmatap.h" 

#include "newmatio.h" 

#ifdef use_namespace 
using namespace NEWMAT; 

#endif 

#define coef 1.0*pow( 10,-6) // define a 


// include.h will get stream functions 
// include.h will get math functions 
// newmatap.h will get include.h 
// need matrix applications 
// need matrix output routines 

// access NEWMAT namespace 


constant coef 


float eigen(float guess_beta, float rap[], float x_node[], float kij[]) 

{ 

//definition of the matrix 
Matrix A( 11,11); 

A<<1<<0<<0<<0<<0<<0<<0<<0<<0<<0<<0 

<<0<<0<<0<<0<<0<<0<<0<<0<<0<<0<<0 

<<0<<0<<0<<0<<0<<0<<0<<0<<0<<0<<0 

<<0<<0<<0<<0<<0<<0<<0<<0<<0<<0<<0 

<<0<<0<<0<<0<<0<<0<<0<<0<<0<<0<<0 

<<0<<0<<0<<0<<0<<0<<0<<0<<0<<0<<0 

<<0<<0<<0<<0<<0<<0<<0<<0<<0<<0<<0 

<<0<<0<<0<<0<<0<<0<<0<<0<<0<<0<<0 

<<0<<0<<0<<0<<0<<0<<0<<0<<0<<0<<0 

<<0<<0<<0<<0<<0<<0<<0<<0<<0<<0<<0 

« 0 << 0 « 0 « 0 « 0 « 0 « 0 « 0 « 0 « 0 « 0 ; 
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A(l,l)=sm(guess_beta*rap[l]); 

A(l,2)=-sin(guess_beta*rap[2]*x_node[l]/x_node[2]); 

A( 1 ,3)=-cos(guess_beta*rap[2] *x_node[ l]/x_node[2]); 

A(2,l)=kij[l]*cos(guess_beta*rap[l]); 

A(2,2)=-cos(guess_beta*rap[2]*x_node[l]/x_node[2]); 

A(2,3)=sin(guess_beta*rap[2]*x_node[l]/x_node[2]); 

A(3,2)=sin(guess_beta*rap[2]); 

A(3,3)=cos(guess_beta*rap[2]); 

A(3,4)=-sin(guess_beta*rap[3]*x_node[2]/x_node[3]); 

A(3,5)=-cos(guess_beta*rap[3]*x_node[2]/x_node[3]); 

A(4,2)=kij[2]*cos(guess_beta*rap[2]); 

A(4,3)=-kij[2]*sin(guess_beta*rap[2]); 

A(4,4)=-cos(guess_beta*rap[3]*x_node[2]/x_node[3]); 

A(4,5)=sin(guess_beta*rap[3]*x_node[2]/x_node[3]); 

A(5 ,4)=sin(guess_beta*rap [3] ) ; 

A(5,5)=cos(guess_beta*rap[3]); 

A(5,6)=-siii(guess_beta*rap[4]*x_node[3]/x_node[4]); 

A(5,7)=-cos(guess_beta*rap[4]*x_node[3]/x_node[4]); 

A(6,4)=kij[3]*cos(guess_beta*rap[3]); 

A(6,5)=-kij[3]*sin(guess_beta*rap[3]); 

A(6,6)=-cos(guess_beta*rap[4]*x_node[3]/x_node[4]); 

A(6,7)=sin(guess_beta*rap[4]*x_node[3]/x_node[4]); 

A(7,6)=sin(guess_beta*rap[4]); 

A(7,7)=cos(guess_beta*rap[4]); 

A(7,8)=-sin(guess_beta*rap[5]*x_node[4]/x_node[5]); 

A(7,9)=-cos(guess_beta*rap[5]*x_node[4]/x_node[5]); 

A(8 ,6)=kij[4] *co s(guess_beta*rap [ 4] ) ; 

A(8,7)=-kij[4]*sin(guess_beta*rap[4]); 

A(8,8)=-cos(guess_beta*rap[5]*x_node[4]/x_node[5]); 

A(8,9)=sin(guess_beta*rap[5 1 *x_node[4]/x_node[ 5 |j; 

A(9,8)=sin(guess_beta*rap[5]); 

A(9,9)=cos(guess_beta*rap[5]); 

A(9,10)=-sin(guess_beta*rap[6]*x_node[5]/x_node[6]); 

A(9,ll)=-cos(guess_beta*rap[6]*x_node[5]/x_node[6]); 

A(l(),8)=kij[5]*cos(guess_beta*rap[5]); 

A(10,9)=4rij[5]*sin(guess_beta*rap[5]); 

A(10,10)=-cos(guess_beta*rap[6]*x_node[5]/x_node[6]); 

A( 10, 1 l)=sin(guess_beta*rap[6]*x_node[5]/x_node[6]); 

A( 1 1 , 1 0)=co s(guess_beta*rap [6] ); 

A( 1 1,1 l)=-sin(guess_beta*rap[6]); 


return A.LogDetenninantO.Value(); //computes the determinant 



void zbrak(float (*fx)(float x, float rap[],float x_node[], float kij[]), float xl, float x2,int n , float 
xbl[], float xb2[], int *nb,float rap[],float x_node[], float kij[]) 

{ 

int nbb=0,i; 
float x,tp,fc,dx; 
x=xl; 

dx=(x2-xl)/n; 

fp=(*fx)(xl ,rap,x_node,kij); 


for (i=l;i<=n;i++) 


fc=(*fx)(x+=dx,rap,x_node,kij); 
if (fp*fc<=0.0) 

{ 

xb 1 [++nbb] =x-dx; 

xb2[nbb]=x; 

if (*nb==nbb)return; 


±TP=fc; 

} 

*nb=nbb; 

} 

//--- 


void find_root (float (*fx)(float x, float rap[], float x_node[], float kij[]), float xbl [J, float xb2[], float 
beta[],int nb, float rap[], float x_node[], float kij[]) 


for (int i=l;i<=nb;i++) 


float xl=xbl[i]; 

float x2=xb2[i]; 

float x_mid, estiinated_root; 

float precision_x = pow( 10,-5); 

int number_iteration=0,max_iteration=50; 

if (fabs((*fx)(xl ,rap,x_node,kij))<=precision_x) 
estimated_root=xl ; 

else if (fabs((*fx)(x2,rap,x_node,kij))<=precision_x) 
estiinated_root=x2; 

else if ((*tx)(xl,rap,x_node,kij)*(*fx)(x2,rap,x_node,kij)>0) 

{ 

cout«"There’s no root between xl and x2.\n"; 
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exit (EXIT_F AILURE) ; 


} 

else if ((*fx)(xl,rap,x_node,kij)*(*fx)(x2,rap,x_node,kij)<0) 

{ 

do 


x_mid=0 . 5 * (x 1 +x2) ; 

if ((*fx) (x_mid,rap,x_node,kij)==0) 
estimated_root=x_mid; 

else if ((*fx)(xl ,rap,x_node,kij)*(*fx)(x_mid,rap,x_node,kij)<0) 
x2=x_mid; 

else 

xl=x_mid; 


number_iteration+= 1 ; 

} 

while (number_iteration<=max_iteration&fabs(0.5*(x2- 
xl))>precision_x&fabs((*fx)(x_mid,rap,x_node,kij))>precision_x); 
estimated_root=x_mid; 

} 

beta[i]=estimated_root; 

} 

} 

// 

float cprime(int i, float x_node[], float K[]) 

{ 

float result=0.0; 

if (1=1) 

return 0; 

else 

{for (int j=2;j<=i;j++) 
result+=x_node[j- l]*(l/K[j- 1]- 1/K[j]); } 
return result; 

I 


// 

float Nm(float A[], float B[], float K[], float x_node[], float alpha[], float b) 

{ 

float sum_l=0.0; 
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for(int i=l;i<=6;i++) 


{ 

float u=2*x_node[i- 1 ] *b/sqrt(alpha[i]); 
float v=2*x_node[i] *b/sqrt(alpha[i]); 

sum_l+=K[i]/alpha[i]*(0.5*(A[i]*A[i]+B[i]*B[i])*(x_node[i]-x_node[i-l])+ 

0.5/b*A[i]*B[i]*sqrt(alpha[i])*(cos(u)-cos(v))+ 

sqrt(alpha[i])/(4.0*b)*(A[i]*A[i]-B[i]*B[i])*(sin(u)-sin(v))); 

} 

return sum_l; 

} 

// 

float Cm(float A[], float B[], float K[], float x_node [], float alpha[], float b, float q) 

{ 

float sum_2=0.0; 
for (int i=l;i<=6;i++) 


float u=x_node[i- 1] *b/sqrt(alpha[i]); 
float v=x_node[i] *b/sqrt(alpha[i]); 

sum_2+=-q/(sqrt(alpha[i])*b*b)*((A[i]*cprime(i,x_node,K)*K[i]*b-B[i]*sqit(alpha[i]))* 
(cos(u)-cos(v))-(A[i]*sqrt(alpha[i])+B[i]*cprime(i,x_node,K)*K[i]*b)* 
(sin(u)-sin(v))+A[i]*b*(x_node[i-l]*cos(u)-x_node[i]*cos(v))- 
B[i] *b*(x_node[i- 1] *sin(u)-x_node[i] *sin(v))); 


return ( sum_2/N n( A, B , K, x_node, alpha, b) ) ; 

} 

// 

// Beginning of the main program 

int main() 

{ 


// Initialization 
float t=0.02;//time (ms) 

float x=39.4*pow(10,-6);//location in the detector 
float q=1.0;//heat flux q=lW/m2.K 

float K[7]={0,71.6,60,0.084,71.6,60,0.209};//layers conductivity (W/m.K) 
float Thickness[7]={0,coef,coef,25.4*coef,coef,coef,10*coef};//layers tluckness (m) 
float Specific_heat[7]={ 0,133. 0,200. 0, 712. 0, 133. 0,200.0, 669. 0};//speciflc heat(J/kg.K) 
float Mass_density[7]={ 0,21450.0, 6880.0, 1289.0,21450.0, 6880.0, 1400.0 };//(kg/m A 3) 
float t_trans=0;//imtialization of the transient temperature 
float phi; 

int layer;// index of layer 
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float alpha[7]={0};//layers diffusivity (m A 2/s) 

float rap[7]={0};//rap gives the ratio x_node[i]/sqrt[alpha[i]] for computing eta[i] 
float x_node[7]={0.0}; 
for (int i=l;i<7;i++) 

{ 

x_node[i]=x_node[i- 1 ]+Tliickness[i] ; 
alpha[i]=K[i]/(Mass_density[i]*Specific_heat[i]); 
rap [i] =x_node [i]/sqrt(alpha[i]); 

} 


float kij [6] = { 0 } ;//kij gives the value of (K[i]/K[i+l])*sqrt(alpha[i+l]/alpha[i]) 

for (i=l;i<6;i++) 

I 

kij[i]=(K[i]/K[i+l])*sqrt(alpha[i+l]/alpha[i]); 

} 

int nb=10; //nb is the maximum of eigenvalues sought 

float beta[l 1]={ 0 } ; 

float xbl[ll];//size always nb+1 

float xb2[ll];//size always nb+1 

zbrak(eigen, 0.0, 290.0, 10, xbl ,xb2,&nb,rap,x_node, kij); 

find root (eigen, xbl, xb2, beta, nb, rap, x_node,kij);//computes the eigenvalues 

//find the layer index corresponding to the value of x 

if (x>0.0 && x<=coef) 
layer=l; 

else if (x>coef && x<=2*coef) 
layer=2; 

else if (x>2*coef && x<=27.4*coef) 
layer=3; 

else if (x>27.4*coef && x<=28.4*coef) 
layer=4; 

else if (x>28.4*coef && x<=29.4*coef) 
layer=5; 

else 

layer=6; 


for (int mdbeta=l;mdbeta<=nb;indbeta++) 

{ 

float b=beta[indbeta]; 
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float A[7]={ 0.0, 1.0}; 
float B[7]={ 0.0, 0.0}; 


for (i=l;i<=5;i++) 


float eta=rap[i]*b; 

float c=x_node [i] *rap[i+ 1] *b/x_node[i+ 1] ; 

A[i+l]=sin(c)*(A[i]*sin(eta)+B[i]*cos(eta))+kij[i]*cos(c)*(A[i]*cos(eta)-B[i]*sin(eta)); 

B|i+l]=cos(c)*(A[i]*sin(eta)+B[i]*cos(eta))-kij[i]*sin(c)*(A[i]*cos(eta)-B[i]*sin(eta)); 


phi=A[layer] * sin(b*x/sqrt(alpha[layer] ))+B [layer] *co s(b*x/sqrt (alpha[layer] )) ; 
t_trans+=exp(-b*b*t)*Cn(A,B,K,x_node,alpha,b,q)*phi;//computes the transient component 

of the temperature rise. 


float t_ss=q*(x/K[layer]+cprime(layer,x_node,K)); //steady- state component. 

//displays on the screen the values of the location, time, and temperature rise 
cout«'\nThe temperature rise at x = "«x«" at t="«t«" is"«(t_trans+t_ss)*pow(10,6) 
«" _micro_K"«"\n\n"; 


return 0; 

} 
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Appendix B: randwalk2d.cpp 


The computer program randwalk2d.cpp determines the temperature at a given location using a two- 
dimensional random walk algorithm. The boundary conditions are specified temperature at the 
upper and lower surfaces, and insulated at the two other surfaces. 


// C++ headers 

#include<io stream. h> 

#include<math. h> 

#include<stdlib.h> 

#include<stdio.h> 

#include<time . h> 

#include<fstream. h> 

// 

List of variables 

N number of random walks 
nx number of nodes in the x direction 

ny number of nodes in the y direction 

L dimension in the x direction, m 

H dimension in the y direction, m 

dx step size in the x direction, m 

dy step size in the y direction, m 

r random number 

(xi,yi) coordinates of the starting point, m 
T_1 temperature of the upper surface, K 
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T_2 temperature of the lower surface, K 


ofstream Point5_file("result.dat");//create a file called result to store the temperature distribution. 

float funct_l (float x, float H) //the upper and lower boundary 
{ return 0*sin(3*x)+0.5; } // of the surface are represented by 

float funct_2 (float x, float H) // trigonometric functions 
{ return 0.0*cos(x)-0. 5; } 


// stait of the mam program 

int main() 

{ 


//Initialization 

clock__t start_time=clock(); 

int N=1500; // number of simulations at each node. 

int nx=100;// number of discrete points in the x-direction 

float L=l().();//dimension in the x direction, m 

float H=l; //dimension in the y direction, y 

int ny=l ()();// number of discrete points in the y-direction 

float dx=L/nx;//step size in the x-direction 

float dy=H/ny;//step size in the x-direction 

float r;//random number; 

float xi=dx/2.0;//x coordinate 

float yi=dy/2.0;//y coordinate 

float T_l=20;//temperature of the upper surface 

float T_2=0;//temperature of the lower surface 

float temp=0.0; 

//change the seed of the random number generator 
int seed=(int) thne(NULL); 
srand(seed); 


for (int i=l;i<=N;i++) 


float x=xi; 
float y=yi; 


while (y<funct_l(x,H) && y>funct_2(x,H)) 


if( x = = ) 
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{ 

x+=dx/2; 

y=y; 

} 


if(x = = L) 


{ 

x-=dx/2; 

y=y; 

} 


r=rand()/(RAND_MAX+1.0); // drawing of the random number 


// search for direction 

if (r<0.25) 


y=y; 

if (x = = dx/2) 
x- = dx/2; 
else 

x- = dx; 


else if(r>=.25 && r<0.5) 

{ 

y = y; 

if (x = =L-dx/2.0) 
x+ = dx/2; 
else 

x+ = dx; 


else if (r>=0.5 && r<0.75) 


{ 

x=x; 

y-=dy; 

I 


else 


x = x; 
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y+=dy; 


// temperature computation after a boundary node is reached 

if (y>0) 
temp+=T_l; 
else 


temp+=T_2; 


//write in the file result 

Point5_file«xi«" "«yi«" "«temp/N«"\u"; 

//compute the machine time 
clock_t stop_time=clock(); 

cout«"Time taken = "«(stop_time-start_time)/CLOCKS_PER_SEC«" secs.Xn"; 
return 0; 

} 
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Appendix C: effcond.cpp 


This code uses a finite difference scheme to determine the temperature profile in a two-dimensional 
geometry. The temperature distribution obtained is used to approximate the heat flux in the y 
direction for the geometry in Figure 19. Equation 3-15 is then used to compute the ratio of the 
effective and bulk conductivities. 


//C++ headers 

#include<io stream. h> 

#include<math. h> 

#include<stdio.h> 

#include<fstream. h> 

//creation of different files 

ofstreamPointl_file("fd2d.dat"); //temperature file 

ofstream Point2_file("x.dat"); //x coordinate 

ofstream Point3_file("y-dat"); //y coordinate 

ofstream Pomt4_fileCh_ratio.dat"); //ratio R/(R+H) and k_eft/k_bulk 


const int ii=5;//number of control volumes in the x_dhection 

const int kk=20;//number of control volume by which R has been divided 


//The function convergence checks if the convergence criteria during the computation of the 
// temperature is satisfied 
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float convergence (float (*M_l)[ii+3], float (*M_2)[ii+3],int jj) 

{ 

float sum=0.0; 

float (*M)[ii+3]=new float [jj+kk+3] [ii+3] ; 
for(int j= 1 ;j<=jj+kk+2;j++) 
for (int i=l;i<=ii+2;i++) 

I 

M[j] [i]=M_2[j] [i]-M_l [j] [i]; 
sum+=M[j][i]; 

} 

return fabs( sum)/( (jj+kk) *ii) ; 

} 

//beginning of the main program 

int main() 

{ 


//Initialization 

int jj=l; //number of control volumes in the y direction 
float (*Tj)[ii+3]=new float [jj+kk+3][ii+3]; 
float L=0.5;//length of the surface 
float R=0.08;// height of the roughness 
float H_ini=0.1; 

float H=H_ini;//initial value of H 

float k_bulk=405.46;//bulk conductivity (W/m.K) 

float x[ii+3]={0); 

float y_step=H_ini/4;//constant step size in the v-direction. 

//the step size being constant in the y-direction, when R increases 

//the number of control volume jj increases also. This keeps the size of the control volumes 
//constant 

for( int step=l;step<=10;step++)// this loop set the number of values of H 
{ //that will be considered, R being constant. 


H+=0.2; 

jj=H/y_step; 

float total_heat=0.0; 
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float (*k)[ii+3]=new float [jj+kk+3][ii+3];//theniial conductivity matrix 


for (int j= 1 ;j<=jj+ 1 ;j++) 

for (int i=l;i<=ii+2;i++) 
k[j] [i]=k_bulk; 


for 0'=jj+2;j<=jj+kk+2;j++) 
for (int i=l;i<=ii+2;i++) 


if (i%2= =0) // specifying this high value for the thermal conductivity 

k[j][i]=pow(10,19); // sets the temperature of those control volumes equals 
else // to the surface temperature and models the geometry in Figure 18. 

k[j][i]=k_bulk; // the remaining control volumes are assigned the bulk conductivity. 


float dx[ii+3|={0}; 

for (int i=2;i<=ii+l;i++) //step size in the x-direction 
dx[i]=L/ii; 


float *dy=new float [jj+kk+3]; //initialization of the step size in the y direction 
dy[0|=0, dy[l]=0, dyQj+kk+3]=0; 


for ( j=2;j<=jj+l;j++) 
dy[j]=H/jj; 


for ( j=]j+2;j<=ij+kk+l;j++) 
dy[j]=R/kk; 


//creation vectors of the coordinates of the control volumes 
//x_coordinates 

for (i=2;i<=ii+2;i++) 

x[ i] =x[i- 1 ] +0. 5 * (dx[i- 1 ] +dx[i] ); 


//y_coordinates 

float (*y)=new float[ jj+kk+3 1; 
for (j= 1 ; j<=jj+kk+2 ; j++) 
y U] =y U- 1 ] +0 • 5 * (dy [j- 1 ] +dy Lj] ) ; 


// 

//Initialization of the matrices 

//These matrices are defined by Vick [1998]. 


effcond.cpp 


93 




float (*aE)[ii+3]=new float[jj+kk+3][ii+3]; 
float (*aW)[ii+3]=new float[ jj+kk+3 1 [ ii+3 1 ; 
float (*aN)[ii+3]=new float [jj+kk+3 1 [ ii+3 1 ; 
float (*aS)[ii+3]=new float [jj+kk+3] [ii+3] ; 
float (*a)[ii+3]=new float [jj+kk+3 1 [ ii+3 1 ; 
float (*b)[ii+3]=new float[jj+kk+3] [ii+3] ; 
float (*p)[ii+3]=new float [jj+kk+3] [ii+3] ; 
float (*q)[ii+3]=new float [jj+kk+3] [ii+3] ; 

// 

//Boundary conditions 

// at x=0 //insulated 

for(j= 1 ; j<=jj+kk+2 ; j++) 


aE [j] [ 1 ] =2*k[j] [2]/dx[2] ; 
a[j] [ 1 ] =aE [j] [ 1 ] ; 
b[j][l]=0; 

} 


//at x=L //insulated 

for(j= 1 ;j<=jj+kk+2;j++) 


a W[j] [ii+2]=2*k[j] [ii+ 1 ]/dx[ii+ 1 ] ; 
a[jj [ii+2] =a W[j] [ii+2] ; 
b[j] [ii+2] =0.0; 

} 


// at y=0 //specified temperature 

for(i= 1 ;i<=ii+2;i++) 

{ 

aN[l][i]=0; 

a[l][i]=l; 

b[l][i]=0.0; //temperature specified is 0 

} 


// at y=H //specified temperature 

for(i= 1 ;i<=ii+2;i++) 

{ 

aS [jj+kk+2] [i]=0; 
a[jj+kk+2] [i] = 1 ; 
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b[jj+kk+2] [i] = 1 00; //specified temperature 1 00 K 

} 

// 

// Computation of the different coefficients aE,aN,aW,aO,a,b 


for (i=2;i<=ii+l;i++) 


for(j=2;j<=jj+kk+ 1 ;j++) 


aW|j][i]=2*k[j] [i-l]*k[j] [i]*dy[j]/(dx[i-l]*k[j][i]+dx[i]*k[j][i-l]); 
aE[j] [i] =2*k[j] [i] *k[j] [i+1 ] *dy[j]/(dx[i] *k[j] [i+ 1 ] +dx[i+ 1 ] *k[j] [i] ); 
aS[j] [i]=2*k[j-l][i]*k|j] [i]*dx[i]/(dy[j-l]*k[j][i]+dy[j]*k[j-l][i]); 
aN[j] [i]=2*k[j] [i] *k[j+ 1] [i] *dx[i|/(dy|j|*k|j+l ||i|+dy[j+l |*k[j] [i]); 
b [j] [i] =0 ;//no source term 
a[j] [i]=aW|j] [i]+aE[j] [i]+aS[j] [i]+aN[j] [i]; 

} 

} 


// 

//initialization of the temperature at all the nodes 
float (*T_guess)[ii+3]=new flo at [jj+kk+3] [ii+3]; 


for(j= 1 ;j<=jj+kk+2;j++) 
for (int i=l;i<=ii+2;i++) 
T_guess[j][i]=0.0; 

float (*T)[ii+3]=new float[jj+kk+3] [ii+3] ; 


for(j= 1 ;j<=jj+kk+2;j++) 
for (int i=l;i<=ii+2;i++) 

Tfj] [i]=T_guess[j] [i] ; 

float (*Ti)[ii+3]=new float[jj+kk+3][ii+3]; 
float (*bi)[ii+3]=new flo at [jj+kk+3 ] [ii+3 ] ; 
float (*bj)[ii+3]=new flo at [jj+kk+3 ] [ii+3 ] ; 


for 0=2;j<=i)+kk+ 1 ;j++) 

for (int i=l;i<=ii+2;i++) 

bi[j] [i] =b[j] [i] +aS [j] [i] *T[j- 1 ] [i] +aN[j] [i] *T[j+ 1 ] [i] ; 
float conv_crit=l; //initialization of the convergence criteria 

// determination of the temperature by iterations in the x and the y direction by using the properties 
// of tridiagonal matrices 
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while (conv_crit>pow(10,-4)) 

{ 

//sweep in the y_direction 


for (j=2;j<=jj+kk+ 1 ;j++) 


P [j] [ 1 ] =aE [j] [ 1 ]/a[j] [ 1 ] ; 
q[j] [l]=b[j] [l]/a[j] [1] ; 


for (i=2;i<=ii+2;i++) 


p[j] [i]=aE[j] [i]/(a[j] [i]-aW[j] [i] *p[j| [i- 1]); 
q[j] [i]=(bi[j] [i]+aW[j] [i] *q[j] [i- 1 ])/(a[j] [i]-aW[j] [i] *p[j] [i- 1]); 
} 


Tfj] [ii+2]=q[j] [ii+2] ; 


for (i=ii+ 1 ;i>= 1 ;i=i- 1 ) 

T|j] [i]=p[j] [i] *T[j] [i+l]+q[j] [i] ; 

for (iiit u=2;u<=jj+kk+l;u++) 
for (i=l;i<=ii+2;i++) 

bi[u] [i]=b[u] [i]+aS[u] [i] *T[u- 1 ] [i]+aN[u] [i] *T[u+l j [i]; 

} 


for(j= 1 ;j<=jj+kk+2;j++) 

{ for(i=l ;i<=ii+2;i++) 

{ Ti[j] [i] =T[j] [i] ; } } 


//sweep in the x_direction 


for (i=2;i<=ii+l;i++) 

forO'= 1 ;j<=jj+kk+2;j++) 

bj[j] [i]=b[j] [i]+aW[j] [i] *T[j] [i- 1 ]+aE[j] [i] *T[j] [i+1] ; 


for (i=2;i<=ii+l;i++) 


p[l][i]=aN[l][i]/a[l][i]; 
q[ 1 ] [i]=b[ 1 ] [i]/ a[ 1 ] [i] ; 


for (j=2;j<=jj+kk+2;j++) 

{ 
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p[j] [i]=aN[j] [i]/ (a[j] [i] -aS [j] [i] *p[j- 1 ] [i]); 

q[j] [i] =(bj[j] [i] +aS [j] [i] *q[j- 1 ] [i] )/ (a[j] [i] -aS [j] [i] *p [j- 1 ] [i] ) ; 


} 

TQj+kk+2 1 [ i |=q | jj+kk+2 1 [i | ; 


f°r(j=jj+kk+ 1 ;j>= 1 ;j=j- 1 ) 

Tfj] [i]=p[j] [i] *T[j+l] [i]+q[j] [i] ; 


for(j=l ;j<=jj+kk+2;j++) 

for (int v=2;v<=ii+l;v++) 

bj[j] [v]=b[j] [v]+aW[j] [v]*T[j] [v-l]+aE[j] [v]*T[j] [v+1]; 


for(j= 1 ;j<=jj+kk+2;j++) 
for(i= 1 ;i<=ii+2;i++) 

Tj[j] [i]=T[j] [i]; 

conv_crit=convergence (Tj,Ti,jj); //compute the convergence value 
} 


Tj [jj+kk+2] [ 1 ] = 1 00; 

Tj Qj+kk+2] [ii+2] = 1 00; 

for(j= 1 ;j<=jj+kk+2;j++) 

I 

for(i= 1 ;i<=ii+2;i++) 

Point l_file«Tj[j] [i] «" 

Point l_file«'\n"; 

} 

for (i=l;i<=ii+2;i++) 

Point2_file«x[i]«"\n";// write the value of x in a file 
for 0= 1 ;j<=jb+kk+2;j++) 

Point3_file«y[j]«"\n";// write the value of y in a file 
for (i=2;i<=ii+l ;i++J 

total_heat+=k_bulk*dx[i]*2/dy[2]*(Tj[2][i]-Tj[l][i]); 

//computation of ratio = k_eff/k_bulk 

float rat io=total_heat*(H+ampl/2)/(L*k_bulk*(b[ jj+kk+2 1 [i]-b[ 1 ] [i])); 
Point4_file«R/(R+H)«" "«ratio«"\n"; // write the values of R/(R+H) and ratio in a file. 
} 

return 0; 

} 
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Appendix D: heterog.cpp 


The program heterog.cpp is adapted from Vick and Scott [1998]. It computes the one-dimensional 
transient temperature distribution of the thermal radiation detector with doped thermal impedance, 
using the two-step process as developed in Section 4-1. 

//Formulation 


T' . - T 


i-1 

m,j 


At 


= K- 


T' . - 2T’ . + T* . 1 

m,j-l m,j m,j+l 


(Ax)^ 



(1) 


rpi rpi— 1 

PJ R ‘ 


— = hTt' -TV 

V rn,j P,J 


p At 
//j control volume. 

//T m j matrix temperature at the j location 
//T P j particle temperature at the j location. 

//The integration of the two equations above yields: 
//From equation (2) we get 

HAt 


C 

rpi P 

Pd - C + HAt 


T i_1 + 


-V . 
C + HAt 

p 


//By replacing T 1 .by (3) the equation (1) becomes: 
P?J 


//&jT m,j aW.T m,j-i aEjT m,j+i - tn 


(2) 


(3) 


(4) 
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HC Ax 

p 


C Ax 

aO -=^7— 
At 


II 

aE = 
j 



aOp,j = 


C +HAt 

p 


aW. = 

j 


k 

W 

Ax 

W 


a. =aO . +aW. +aE. +aO b.=aO .sT 1 f + aO .*T J 1 

J m-j j j P,J J m,J m,j p,j p,j 


// C++ headers 

#include<math . h> 

#include<iostream . h> 

#include<f stream . h> 

#include<stdio . h> 

#define coef pow(10,-6) 

//creation of data files 

ofstream Point l_file ( "Matrix_temp . dat ") ; //matrix temperature file 
of stream Point2_f ile ( "Particles_temp . dat" ) / // particles temperature file 
ofstream Point3_f ile ( "y . dat" ) ; //location 

ofstream Point4_fileCTemp_act_junct.dat"); //active junction temperature file 

const int t_cvy=60; //total cv in the y_direction=sum of ny 

//main program 

int main ( ) 

{ 

//initialization 

double density=0 . 01; //volume fraction of the particles 
double dp=pow ( 10 , -7 ) ; //diameter of the particles, m 
double rocp=1860 . 2*pow (10, 3) / //heat capacity the particles 
double vp=3 . 14*pow (dp, 3) /6 . 0/ //volume of a particle 

double Np=density/vp; //number of particles in the thermal impedance layer per 
volume; 

double h=10; contact conductance, W/nr.K 

double Ap=3 . 14*dp*dp; //surface area of the particles, irT 

float Q=pow ( 10 , 6) ; //heat flux, W/m 2 

double tmax=0.1;// time maxi, s 

int n time=100; // number of time steps 

int count=l; 

double dt=tmax/n time;// time step, s 
double t=0; //time 

double kp=l . 5 9; //thermal conductivity of the particles, W/m,K 

double H [7] ={ 0, coef , coef , 25 . 4*coef, coef , coef , 10*coef} ; //thickness of the layers 
double conduct [ 7 ] ={ 0 , 71 . 6, 60 , 0 . 084 , 71 . 6, 60 , 0 . 209 } ; // layers conductivity 
(W/m.K) 

double ny [7 ] ={ 0 , 10 , 10 , 10 , 10 , 10 , 10 } ; //number of control volumes of each layer 
double np [ 7 ]={ 0 , 0 , 0 , Np, 0 , 0 , 0 }; //number of particles 

double Specific heat[7]={0, 133. 0,200. 0,712. 0,133. 0,200. 0,669.0};// specific 
heat ( J/kg . K) 

double Mass density [7 ] = { 0 , 21450. 0,6880. 0,1289. 0,21450. 0,6880. 0,1400.0};/ /Mass 
density of the layers (kg/m A 3) 
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int j j [ 7 ] = { 0 } ; 

double bio=h*dp/ ( 6*kp) ; //computes the Biot number 
double roc [ 7 ] = { 0 } ; //heat capacity of the layers 
for (int i=l; i<=6; i++) 

{ 

roc [i]=Mass density [i] *Specific heat[i]; 
j j [i] = j j [ i— 1 ] +ny [i] ; 

} 

double km[t cvy+3] = { 0 . 0 } ; 
double Rocm[t cvy+3] ={ 0 . 0 } ; 
double K[t cvy+3 ]={ 0 . 0 } ; 
double Cm[t cvy+3] = { 0 . 0 } ; 
double Cp[t cvy+3] ={ 0 . 0 } ; 

int c=l; 

for (int m=l/m<=6/m++) 

{ 

for (int j=c; j<= ( j j [m] +2 ) ; j++) //effectives properties of the layers of the 

//thermal detector 

{ 

km [ j ] =conduct [m] ; 

Rocm [ j ] =roc [m ] ; 

K [ j ] = ( 1-np [m] *vp) *km[ j] / 

Cm [ j ] = (1-np [m] *vp) *Rocm [ j ] ; 

Cp [ j ] =np [m] *vp*rocp; 

} 

c=j j [m] +2/ 

} 

double dy [t_cvy+3] = { 0 . 0 } / //step size vector, m 
double y [t_cvy+3 ] = { 0 . 0 } ; //location vector, m 

int d=2 ; 

for (m=l;m<=6;m++) 

{ 

for (int j=d; j<= j j [m] +1; j++) 
dy [ j]=H[m] /ny [m] / 
d=j j [m] +2; 

} 

for (int j=2;j<=t cvy+2;j++) 

Y [ j]=y [ j-1] +0 • 5* (dy [ j-l] +dy [ j] ) / 

for ( j=l; j<=t_cvy+2; j++) 

Point3 f ile«y [ j ] <<" \n" ; 

//Initial temperature 

double Tmprev[t cvy+3] ={ 0 . 0 } ; //matrix 
double Tpprev[t cvy+3] ={ 0 . 0 } ; //particles; 

Point4 file«0«" "«0«" \n" / //initial value of the active junction temperature. 

double aN[t cvy+3] ={ 0 . 0 } / 
double aS [t cvy+3 ]={ 0 . 0 } / 
double a[t cvy+3 ] = { 0 . 0 } ; 
double b[t cvy+3 ] ={0.0}/ 


hetero.cpp 


100 



double p[t cvy+3 ] = { 0 . 0 } ; 
double q[t cvy+3 ]={ 0 . 0 } ; 
double amO[t cvy+3] ={ 0 . 0 } ; 
double apO[t cvy+3] ={ 0 . 0 } ; 
double Tm[t cvy+3] ={ 0 . 0 } ; 
double Tp[t cvy+3 ] ={0.0}; 


//Boundary conditions 

// @ y=0 specified temperature, 0 K 

aN [1] =0; 
a [ 1 ] = 1 ; 
b [ 1 ] = 0 ; 

// @ y=H specified heat flux Q 

aS [t cvy+2 ] =2*km [t cvy+l]/dy[t cvy+1]; 
a [t_cvy+2 ] =aS [t_cvy+2 ] ; 
b [t_cvy+2 ] =Q; 


//Computation of the different coefficients aN, aS, amO, apOa, b 

for ( j=2 ; j<=t_cvy+l ; j++) 

{ 

aS [ j ] =2 *K [ j — 1 ] *K[j] / (dy [ j -1 ] *K [ j ] +dy [ j ] *K [ j-1 ] ) ; 
aN [ j ] =2*K [ j ] *K [ j + 1 ] / (dy [ j ] *K [ j + 1 ] +dy [ j + 1 ] *K [ j ] ) ; 
amO [ j ] =Cm [ j ] *dy [ j ] /dt ; 

apO [ j ] =Cp [ j ] *h*Np*Ap*dy [ j ] / (Cp [ j ] +h*Np*Ap*dt) ; 
a [ j ] =amO [ j ] +apO [ j ] +aS [ j ] +aN [ j ] ; 

} 


//computation of the matrix and particle temperature by solving the coupled 
//Equations 3 and 4 . Equation 4 is solved by using the tridiagonal matrix 
properties . 

for (count=l; count<=n_time; count++) 

{ 

t+=dt ; 


for ( j=2 ; j<=t_cvy+l ; j++) 

b [ j ] =amO [ j ] *Tmprev [ j ] +apO [ j ] *Tpprev [ j ] ; 

P [ 1 ] =aN [ 1 ] /a [ 1 ] ; 

q [ 1 ] =b [ 1 ] / a [ 1 ] ; 


for (j=2;j<=t cvy+2 ;j++) 

{ 

p [ j] =aN [ j] / (a [ j ] -aS [ j] *p [ j-1] ) ; 

q[ j]=(b[ j]+aS[ j]*q[ j-1] ) / (a[ j]-aS[ j]*p[ j-1] ) ; 


Tm[t cvy+2 ]=q[t cvy+2]; 

for ( j=t_cvy+l ; j>=l; j= j-1) 

Tm [ j ] =p [ j ] * Tm [ j + 1 ] +q [ j ] ; 
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for (j=l;j<=t cvy+2;j++) 

{ 

Tmprev [ j ] =Tm [ j ] ; 

Pointl file«t*pow (10, 3) «" "<<Tm [ j ] «" \n" ; 

} 

//write the temperature of the active junction in a file 
Point4 file«t*pow (10, 3) «" "«Tm [ j j [4] +1] «"\n" ; 


for ( j=j j [2] +2; j<=j j [3] +1; j++) 

{ 

Tp [ j ] = (h*Np*Ap*dt*Tm [ j ] +Cp [ j ] *Tpprev [ j ] ) / (Cp [ j ] +h*Np*Ap*dt ) ; 
//write the temperature of the particles in a file 
Point2_f ile«t«" "«Tp [ j]«"\n"; 

} 


Point5 file<<t*pow (10, 3) <<" "«Tm [ j j [2 ] +4 ] «" "«Tp [ j j [2 ] +4 ] «" \n" ; 

for ( j=l; j<=t_cvy+2; j++) 

Tpprev [ j ] =Tp [ j ] / 


} 

return 0/ 

} 
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